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CHAPTER  I 


INTRODUCTION 


Chisholm's  Third  Law: 

Proposals,  as  understood  by  the  proposer,  will 
be  judged  otherwise  by  others. 

Lilly's  Metalaw: 

All  laws  are  simulations  of  reality. 

Levy's  Ninth  Law: 

Only  God  can  make  a  random  selection. 


Both  discrete-event  and  continuous  simulations 
are  frequently  used  to  model  analytically  intractable 
stochastic  systems.  Examples  of  such  systems  are  city 
traffic  flow,  telecommunications  networks,  military  war 
games,  and  computer  time-sharing  systems.  Due  to  the 
complexity  of  these  systems,  simulation  is  often  the 
only  viable  means  of  studying  them.  There  are,  however, 
significant  tactical  problems  involving  autocorrelated 
outputs  and  initialization  bias  which  complicate  the 
analysis  of  such  simulations.  The  regenerative  method 
is  a  recently  developed  technique  for  analyzing  simula¬ 
tion  output  which  promises  to  alleviate  these  problems* 
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1.1  Estimating  Mean  Response  Time  in  the  M/M/1  Queue 
To  illustrate  the  use  of  the  regenerative 
method,  we  will  consider  the  estimation  of  the  expected 
stationary  waiting  time  in  an  M/M/1  queue.  The  pro- 
#  cess  to  be  analyzed  is  {Wj  :  j  >_  l},  the  waiting  times 

of  successive  customers.  In  the  regenerative  method, 
times  at  which  an  arriving  customer  finds  the  system 
empty  and  idle  are  of  particular  interest.  Given  an 
arrival  rate  X,  a  service  rate  p,  and  a  traffic  inten¬ 
sity  p  »  A/p  <  1,  it  can  be  shown  [CRAN74(I)]  that  with 
probability  1  there  is  an  increasing  sequence  of 
integer-valued  random  variables  Bj.,  &2,  B3,  ...  such 
that  Wgj  «  0  for  all  positive  integers  j.  That  is,  any 
customer  numbered  Bj  will  arrive  to  find  the  system 
empty  and  will  be  serviced  immediately.  The  arrival 
times  for  the  customers  indexed  by  Bj  define  regenera¬ 
tion  points  for  the  system.  At  each  of  these  points, 
system  operation  begins  anew  with  the  same  probabilistic 
4  structure  that  governed  it  at  the  time  of  the  first 

arrival.  Thus,  the  process  Wj  may  be  partitioned  into 
blocks  or  cycles  iw^  *  Bj  £  i  <  6j+i>  which  are  indepen¬ 
dent  and  identically  distributed  (iid).  Any  measure¬ 
ments  taken  on  these  cycles  will  also  be  iid.  Let 
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Si+1-1 

Yi  -  E  Wj  and  ai  ■  6^+1  -  @i» 

where  Y*  represents  the  total  waiting  time  of  customers 
served  in  the  ith  block  and  denotes  the  number  of 
customers  in  the  ith  block.  Then  the  random  vectors 
(Y*,  a±)  are  iid.  From  regenerative  analysis,  we  have 


[CRAN75U)] 


W  and  yw  -  E[W]  -  EIYxJ/ElaiJ 


Given  the  observations  (Yj,,  ai)  s  1  <  i  <  n  ofn  simu¬ 
lated  cycles,  we  compute  the  sample  means  Y  and  a  in 
order  to  form  the  following  ratio  estimator  of  u^: 


r  -  Y/a 

By  computing  the  sample  variance  S2  of  the  quantities 
(Yi  -  ?cii)  i  1  <,  i  <,  n  ,  we  also  obtain  an  approximate 
100(1  -  6)%  confidence  interval  for  Uw(CRAN77] * 


r  '  *l-$/2 


If  we  require  that  the  half-length  of  this  interval  fall 
within  some  proportion  d  of  the  point  being  estimated, 
then  the  required  number  of  cycles  is  given  by  [DALS68] 

n  -  (t1.5/2/d)2(p3  -  4p2  +  5p  +  2)/[p(l  -  p)], 

and  the  expected  amount  of  simulated  time  to  obtain  n 
cycles  is 

t  “  Jt/Vd)2(P3  -  4P2  ♦  5P  +  2)/(  P  (1  -  P)2) . 

1-5/2 

Table  1.1  shows  the  required  number  of  cycles  for 
various  traffic  intensities  p  with  A  ■  1  when  we  seek  a 
95%  confidence  interval  whose  half-length  is  5%  of  the 
estimand  ity*  Although  the  M/M/1  queue  is  not  typical  of 
most  real-world  systems  because  it  is  analytically 
tractable,  it  is  still  representative  of  the  com¬ 
putational  costs  required  by  the  regenerative  method  of 
simulation. 

1.2  Problem  statement 

Given  a  regenerative  queueing  simulation,  the 
problem  we  face  is  to  reduce  the  sampling  coats  required 
to  obtain  acceptable  precision  in  simulation-baaed  esti¬ 
mators.  To  do  so,  we  must  choose  applicable  vagiaama 
reduction  techniques  (vxts)  which  provide  increased  pre¬ 
cision  while  maintaining  or  rcduoirg  the  simulation  run 
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TABLE  1.1 

REQUIRED  SAMPLING  VOLUME  FOR  M/M/1  WAITING  TIMES 
WITH  «  -  .05,  d  -  .05,  X-  1.0 


Praffic 


Number  ol 


Expectec 


Intensity 

p 

Cycles 

n 

Run  Length 
t 

.01 

318,132 

321,345 

.05 

72,469 

76,283 

.10 

42,019 

46,688 

.15 

32,099 

37,764 

.20 

27,353 

34,191 

.25 

24,714 

32,952 

.30 

23,174 

33,106 

.40 

21,923 

36,538 

.50 

22,282 

44,563 

.60 

24,178 

60,445 

.70 

28,413 

94,711 

.80 

37,955 

189,776 

.90 

68,107 

681,073 

.95 

129,316 

2,586,327 

.99 

620,849 

62,084,898 

length.  In  applying  these  VRTs,  we  are  constrained  by 
the  requirement  for  valid  confidence  intervals.  Much 
recent  research  has  been  devoted  to  VRTs  in  this 
setting.  Most  of  this  research  has  produced  mixed 
results. 

1.3  Objectives  of  the  Research 

The  goals  of  this  research  are  the  development , 

.  implementation,  and  evaluation  of  techniques  to  improve 
the  efficiency  of  regenerative  queueing  simulations.  He 
have  chosen  to  restrict  ourselves  to  VRTs  which  do  not 
alter  the  sample  path  generated  by  a  simulation  model. 

We  sought  to  employ  robust  procedures  which  make  effec¬ 
tive  use  of  auxiliary  information  produced  by  the  simu¬ 
lation  and  which  may  be  applied  to  a  variety  of  real- 
world  queueing  situations.  In  light  of  this,  we  have 
developed  a  regression-based  technique  using  internal 
control  variables  to  achieve  increased  efficiency  in  the 
simulation  of  regenerative  queueing  networks  and  we  have 
established  the  asymptotic  properties  of  the  procedure 
which  ensure  its  validity  and  efficiency.  This  method 
is  evaluated  in  several  selected  experimental  systems  to 
demonstrate  its  potential  value. 


CHAPTER  II 


LITERATURE  REVIEW 

2.1  The  Regenerative  Method 

Crane  and  Iglehart  pioneered  the  regenerative 
method  in  a  series  of  recent  papers  [CRAN74(I),  (II), 
75(111),  (IV),  IGLE75,  76],  The  first  of  these  deals 
with  the  general  multiserver  queue. 

Consider  a  GI/G/1  queue  with  service  rate  y, 
arrival  rate  X,  and  traffic  intensity  p  *  X/y  <  1.  For 
simplicity  we  assume  that  the  zeroth  customer  arrives  at 
time  t  «  0  and  finds  the  system  empty  and  idle.  In  gen¬ 
eral,  the  ith  customer  arrives  at  time  t*,  experiences 
a  waiting  time  Wj_,  and  finally  receives  a  service  of 
duration  v^.  The  interarrival  times  are  given  by 
Uj,  *  ti  -  ti-i,  i  £  1.  We  assume  that  the  two  sequences 
(ui  :  i  >,  1}  and  {vi  :  i  _>  0}  are  mutually  independent 
and  consist  of  iid  random  variables. 

For  this  system,  variates  of  interest  include 
Q(t),  the  number  of  customers  in  the  system  at  time  t; 
Wi,  the  waiting  time  of  the  1th  customer;  W(t),  the 
workload  the  server  sees  at  time  t;  B(t),  the  amount  of 
time  during  the  period  [0,t]  in  which  the  server  is 
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busy;  and  D(t),  the  number  o£  customers  who  completed 

service  during  the  period  [0,tl.  Let  Xn  *  vn_j_  -  un  and 
n 

Sn  «  So  *  0.  Then  for  the  process 

{Wn  :  n  >_  0},  we  have  the  following  recursive 
relationship: 

W0  -  0  and  Wn  «  (W,^  +  Xn]+,  n  >  1  (2.1) 

where  [a]+  =  max  {0,a}.  It  can  be  shown  that  with  prob¬ 
ability  one  there  exists  a  strictly  increasing  sequence 
of  integer-valued  random  variables  {Bk  :  k  >.  0}  such 
that  Wgk  *  0  for  k  ^  0.  Customers  numbered  {3^}  are 
those  arriving  to  find  the  system  empty  and  idle.  The 
arrival  of  the  customer  with  index  3k  initiates  the  kth 
busy  period  during  which  the  server  remains  occupied. 

At  the  end  of  this  busy  period,  the  kth  idle  period 
begins  and  lasts  until  the  customer  indexed  by  3)^ 
arrives.  A  busy  period  and  its  succeeding  idle  period 
form  a  busy  cycle  or  tour.  The  number  of  customers 
served  in  the  kth  busy  period  is  given  by 
°i  *  3i  -  Sj-i,  i  >  1. 

Using  Crane  and  Iglehart's  notation,  we  let 
*k  •  (vk-l,uic)  and  Yk  *  (ak/^.j+i,  XBk-j+2*  ...»  Xpk). 
The  {Vfc  :  k  £  1}  are  iid  and  allow  the  observations  to 


be  divided  into  iid  cycles.  On  each  cycle  we  will 
observe  the  following: 


3k  -  1 

Hk  *  2  Vi  ,  (2.2) 

i»3k-l 

Bk 

Cjt  *  £  ui  ,  (2.3) 

i-Sk-l+1 

vk  *  Sk  “  nk  »  <2*4> 

<*k-l 

*k(i)  -  £  (Se)c-1+j  -  Sg^)  ,  (2.5) 

j»l 

yk(2)  -  Yk(1)  +  nje  ,  (2.6) 

and 

<*k-l 

Y*(3)  -  e  ns0k_1+j  -  S6k_1)V0k.1+j  + 

j-o 

2  <v6k-l+j)2]  '  (2‘7) 


where  n*,  5k/  and  v*,  respectively,  represent  the 
durations  of  the  busy  period,  busy  cycle,  and  idle 
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period  of  the  kth  regenerative  cycle;  Yfc^2*  and  Yfc*3*, 
respectively,  denote  the  integrals  of  Q(t)  and  W(t)  over 
the  kth  cycle.  Mote  that  all  of  these  random  variables 
form  iid  sequences. 

The  regenerative  property  of  the  Gl/G/1  queue 
insures  that  the  processes  {Wj},  {Q(t)},  and  {W(t)}  con¬ 
verge  in  distribution  to  corresponding  steady-state 
variates: 

wn  ♦  W,  Q( t )  -*•  Q  ,  and  W(t)  ♦  W*  . 

From  stochastic  processes  theory,  we  have: 


*  CO 

Etefc]  -  exp(  E  p  {Sn  >  0}  /n)  , 
n*l 

(2.8) 

Etnx]  -  E(akl/y  , 

(2.9) 

EUkl  "  Ef-jcJ/l  , 

(2.10) 

Elvk]  -  (1  -  p  )E [ajjJ/X  , 

(2.11) 

E[YkU))  -  E[W]E[ak]  , 

(2.12) 

EtYk(2)]  «  (XE  [W]  ♦pJEUid  -  EtQlEUjcl 

,  (2.13) 

and 


££§§Sg&| 
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EtYfcU)]  -  <pE[W]  +  ^EIv£]  )E[5kl 

-  ElW*]E[5ic]  .  (2.14) 

Thus  £rom  our  observed  random  variables,  we  may 
calculate  point  estimates  £or  the  expected  values  of  W, 
Q,  and  W*. 

Crane  and  Iglehart  [CRAN75(lli)]  generalized 
their  analysis  o£  multi-server  queues  to  discrete-event 
simulation  models.  We  now  let  X(t)  denote  the  model 
status  vector  at  (simulated)  time  t.  If  (X(t)  t  t  >  0} 

•v  * ” 

is  a  regenerative  process,  then,  subject  to  mild  regu¬ 
larity  conditions.  Crane  and  Iglehart  showed  that  a 
steady-state  system  status  vector  X  exists  such  that 
X(t)  *  X. 

*v 

If  the  goal  of  the  simulation  is  to  estimate 
some  steady-state  performance  measure  r  *  E[f(X)],  then 
successive  regeneration  epochs  (Sj)  of  the  simulation 
model  are  observed  in  order  to  accumulate  the  following 
measurements  over  each  cycle: 

•  4 

®j+l 

Yj  -  /  f (X( t) )dt  and  -  6 j  . 

•i 


(2.15) 
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The  random  vectors  {(Yj,aj)  :  j  >  1}  are  therefore  iid. 
Crane  and  Iglehart  show  that  under  certain  mild 


regularity  conditions, 


r  -  E  [YjJ  /E  ( <*],] 


(2.16) 


For  some  systems  these  regularity  conditions  are  that 
the  distribution  function  of  is  not  arithmetic  and 
that  0  <  E[axl  <  00 . 

Suppose  that  n  regenerative  cycles  {(Yj,otj)  j 
1  £  j  <  n)  have  been  observed.  Let  Vj  ■  Yj  -  raj.  The 
Vj's  are  iid  with  mean  zero.  If  V,  Y,  and  a  are  the 
sample  means  of  V,  Y,  and  a,  respectively,  then  V  »  Y  - 
ro.  If  0  <  a2  *  Var[Vj  <•,  the  Central  Limit  Theorem 


gives  us 


lim  PJ -  <  z  »  #(z)  for  all  z 

n-H*>  fa  /  /n  ~ 


(2.17) 


where  *  is  the  standard  normal  distribution  function. 
This  result  may  be  written  as 


lim  ,  <  z 

n-H»  c  /  / n 


*(z)  , 


(2.18) 


lim  P  \ - SL “-r_  1<  z  .  $(z)  ,  (2.19) 

n-*K»  )a/(/n  *  a)J 


where  r  »  Y/a.  in  order  to  obtain  a  confidence  interval 
from  (2.19),  we  must  estimate  a.  Now, 

02  -  ElVil  *  EUYx  -  rax)2) 

-  VarlYxl  -  2rCov(Yx,ax]  +  r2var(ax]  .  (2.20) 
2  2 

Let  S‘  and  S‘  be  the  sample  variances  of  Y  and  a, 
respectively,  and  let  S2Q  be  the  sample  covariance  of  Y 
and  a: 

Sy  -  U/(n  -  1)]  E  (Y-j  -  Y)2, 

j-1  J 

2  n 

S‘  -  ll/(n  -  1)]  E  (04  «  a)2  , 

a  j-1 


and 


n  _ 

[l/(n  -  1)]  E  (Yj  -  Y)  (cm  -  a) 2 
j-1  J 


(2.21) 


Let  S2  -  S2  -  2rSya  +  r2S2;  then  we  have  S2  n-«^  cr2  with 
probability  1,  and  it  follows  that 

a  | 

li»  p  — S-Z -g-l<  *  -  ♦(*)  for  all  z  . 

n-Mo  S/(/n  a)  | 


(2.22) 


An  approximate  100(1  -  6)  percent  confidence  interval  is 
then  given  by 

r  +  z,  ./2  *  — (2.23) 
16/2  /n  a 

where  z ^  *  *_1(1  -  6/2). 

Crane  and  Iglehart  [CRAN75( III ) ]  point  out  that 
it  may  be  possible  to  find  a  second  sequence  of  random 
variables  :  i  >_  0>  which  also  define  regeneration 
points.  If  confidence  intervals  with  lengths  I(t)  and 
I'(t),  respectively,  are  constructed  from  the  sequences 
{ (Yi,ai ) }  and  {(Yj.,ai)}  each  of  which  are  based  on  a 
simulation  run  length  t,  then  I(t)/I'(t)  +  1  with 

probability  1.  Thus,  if  a  simulator  may  choose  between 
two  or  more  regeneration  sequences,  the  lengths  of  the 
confidence  intervals  will  be  approximately  the  same  for 
large  t. 

The  preceding  discussion  was  based  upon  the 
assumption  that  a  return  state  can  be  found  with  the 
property  that  the  expected  time  between  returns  is 
finite  and  small  enough  so  that  a  reasonable  number  of 
cycles  will  be  observed  during  the  simulation.  This  may 
not  always  be  the  case.  Crane  and  Iglehart  (CRAM75(IV)1 
offer  some  approximation  techniques  for  obtaining 


confidence  intervals  when  the  simulation  does  not 
contain  a  renewal  process. 

The  first  method  presented  to  deal  with  this 
problem  is  complete  state-space  discretization.  For  a 
queueing  system,  this  technique  requires  the  approxima¬ 
tion  of  the  interarrival  and  service  time  distributions 
which  may  only  take  on  values  {0,  x,  2t,  ...},  x  >  o. 

The  choice  of  x  is  critical  to  this  method.  The  smaller 
x  is  chosen,  the  closer  the  new  process  will  mimic  the 
original.  A  smaller  t  will  also,  however,  result  in 
fewer  observed  cycles  within  a  fixed  simulation  run. 

A  second  method  for  handling  the  regeneration 
problem  is  partial  state-space  discretization.  Consider 
the  customer  waiting  times  {Wn  :  n  >  ll.  We  shall 
modify  the  original  process  as  follows: 

,  d,  if  d  -  e  <  Wn  <  d  +  e 
W_  -  “  (2.24) 

Wn,  otherwise 

where  d  is  a  fixed  waiting  time  and  e  >  0  is  the  half- 
length  of  the  "trapping  interval”  around  d.  The  entry 
times  {et(d>  :  i  >  1}  to  d  are  the  regeneration  times 
for  the  modified  process.  The  confidence  interval 
methods  developed  earlier  may  now  be  applied  to  the 
modified  process.  As  e  0,  the  modified  process 


should  approximate  the  original  process  more  closely, 
but  fewer  of  the  required  regeneration  cycles  will  be 
observed. 

The  third  technique  suggested  by  Crane  and 
Iglehart  is  stochastic  bounding.  This  method  requires 
the  simulator  to  define  two  new  processes  which  bound 
the  original.  Thus,  in  finding  confidence  intervals  for 
the  new  processes,  a  confidence  interval  for  the 
original  may  be  found.  Stochastic  bounding  uses  the 
same  scheme  as  partial  state-space  discretization  in 
that  it  is  baaed  upon  selecting  a  trapping  interval 
[d  -  e,  d  +  ej  about  d.  For  the  customer  waiting  time, 
our  new  processes  are 

,  d-e,  if  d  -  e  <  »„  <  d  +  e 
W  -  ”  (2.25) 

Wn,  otherwise 

and 

.  d+e,  ifd-e<Wn<d+e 
W  -  “  (2.26) 

Wn,  otherwise  . 

These  processes  bound  the  originals 
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for  all  w  >_  0;  and  if  f  :  R+  -*•  R  is  a  measurable, 
monotonically  increasing  function, 

E{f(W')l  <  E[f(W)]  <  E[f(W")]  (2.28) 

V  Jt 

where  W*  ^  W  and  W“  —*■  W".  Using  d  -  e  and  d  +  e  as 
n  n 

the  regenerative  states  for  (W'n}  and  {W^  },  respec¬ 
tively,  we  find  their  100(1  -  6/2)%  confidence  intervals 
in  the  standard  manner.  If  we  take  the  lower  limit  of 
the  confidence  interval  for  E [f (W ) ]  and  the  upper  limit 
of  the  confidence  interval  for  Elf(W")J  to  construct  a 
confidence  interval  for  E[f(W)],  we  have  a  minimum  of 
100(1  -  6)%  coverage. 

The  last  method  suggested  by  Crane  and  Iglehart 
for  approximating  a  regenerative  system  is  the  use  of 
approximate  regeneration  times.  As  in  partial  state 
space  discretization,  the  entry  times  {(3j(d)  s  i  >  1}  to 
a  trapping  interval  Id  - e  ,  d  +  e]  are  taken  as  the 
regeneration  times.  Observations  of  the  process  are 
collected  as  for  any  regenerative  system,  and  confidence 
intervals  are  constructed  by  the  standard  method.  For  a 
small  e,  this  method  produces  results  similar  to  those 
which  would  be  obtained  from  the  original  process  in  the 
normal  manner.  The  blocks  for  this  modified  process  are 
not  iid.  Zf  the  correlation  between  successive  blocks 


is  calculated,  the  simulator  can  get  some  idea  of  the 

validity  of  his  confidence  interval. 

Gunther  and  Wolff  [GUNT80]  recently  proposed  the 

Almost  Regenerative  Method  for  handling  regenerative 

processes  which  have  infrequent  regeneration  points. 

For  simplicity  they  assumed  that  at  each  regeneration 

point  0n  there  is  a  change  of  state  in  the  process 

{ X( t ) }  from  some  fixed  state  u  to  some  fixed  state  v,  u, 

v  e  E,  the  state  space  of  X.  (They  indicate  that  this 

is  not  necessarily  a  requirement,  but  that  it  gives  them 

a  method  to  compare  their  technique  to  the  regenerative 

method. )  Let  U  and  V  be  disjoint  subsets  of  E.  (U  and 

V  do  not  have  to  partition  E. )  Let  0  *  denote  the  time 

n 

of  the  nth  transition  of  (X(t)}  from  any  state  in  0  to 
any  state  in  V.  For  the  Almost  Regenerative  Method, 
the  {0^  :  n  >_  0}  are  the  regeneration  points  and  the 
duration  of  the  nth  cycle  iso'«8'-B'1,n>l.  The 
process  (X(t)}  is  observed  for  m*  "cycles”  and  the  pairs 
<(*{,«!)>  are  collected.  The  estimate  of  B[f(X)]  is 
given  by 

r*  -  Y'/S“*  .  (2.29) 

If  the  underlying  process  (X(t)>  is  regenerative  and 
u  e  0  and  v  e  V,  Gunther  and  Wolff  were  able  to  show 


that  r1  is  a  consistent  estimator  of  E[f(X)]  with  an 
asymptotically  normal  distribution. 

Gunther  and  Wolff  offered  total  amount  of  work 
in  the  network  and  number  in  system  as  methods  of 
selecting  U  and  V.  Their  selection  is  based  upon 
desiring  to  estimate  specific  response  variables.  The 
choice  of  an  appropriate  U  and  V  is  critical,  and  an 
improper  selection  may  give  uncertain  results.  For 
well-chosen  U  and  V,  they  were  able  to  obtain  smaller 
variances  than  the  standard  regenerative  method, 
although  they  point  out  that  this  may  be  directly 
attributed  to  the  larger  number  of  cycles  which  they 
were  able  to  obtain. 

2.2  Variance  Reduction  Techniques 

We  have  seen  that  large  sample  sizes  are 
frequently  required  in  regenerative  simulation.  We 
shall  next  examine  some  common  variance  reduction 
techniques  which  may  be  employed  to  reduce  run  length. 

2.2.1  Stratified  Sampling 

In  stratified  sampling,  observations  are 
collected  on  the  simulation  response  Y  and  a  stratifica¬ 
tion  variable  X  which  has  a  known  distribution  and  a 
strong  but  highly  nonlinear  dependence  on  Y.  Strata 


•  m 
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* 

.  n 


are  formed  by  partitioning  the  range  of  X  Into  L  strata 
{Sfc  :  1  <  k  <  L).  Bach  stratum  ia  then  examined 
separately  to  obtain  a  corresponding  stratified  random 
sample  of  the  response  Y.  For  1  <  Is  <  L  let 


W*  «  P{XeSk}  , 


Uyk  *  Elfj  XeSfc]  # 


(2.30) 

(2.31) 


and 


o^y.  *  E  ( (Y  -  wyic)2|  xesjcl 


(2.32) 


Suppose  we  sample  n*  simulation  responses  {  Y*j  * 

1  <  j  <  1%}  from  each  stratum  S*.  The  sample  mean  for 
each  stratum  is  given  by 


nk 

*k  -  U/nk)  £  *kj 
j*l 


and  the  stratified  estimator  of  yy  is 


ft 


£  ltfcffc 


f.  ia  am  unbiased  estimator  and  ham  variance 


Var[Ys]  *  E  wv°vv/nfc  • 
k-1  K  x* 


(2.35) 


I£  sampling  over  the  strata  is  done 
proportionally,  i.e.,  n*  »  nWjc,  1  <,  k  <  L,  then  Y3p 
denotes  the  stratified  estimator  o£  yy;  and  we  have: 


Var(YSp]  -  (l/n^E^oy* 


2  u  7 

Oy/n  -  (1/n)  E  W^lyy^  -  y y ) 
k-1 


(2.36) 


If  the  stratum  means  are  not  all  equal,  Var[Yap]  will  be 

2 

strictly  less  than  the  variance  Uy/n  of  the  mean  of  the 

unstratified  sample.  If  the  within-stratum  variances 
2 

{Cy^}  widely  differ,  it  is  desirable  to  sample  more 
heavily  from  the  strata  with  larger  variability.  Using 
the  optimal  (Neyman)  sampling  scheme 


nk  -  n  *  (WkayfcJ/t  E  W.Oy,)  ,  1  <  k  <  L  ,  (2.3' 

£«1 


the  maximum  variance  reduction  ia  achieved.  This 
allocation  presupposes  knowledge  of  the  {Oy^}.  In  most 
simulations  this  is  unknown  but  may  be  estimated  in 


£*0*0 «  pK  <90 


23 

preliminary  runs.  Once  the  allocations  have  been 
established,  the  sample  variances 

nk 

*  [l/(njt  -  1)3  I^OTkj  -  Yk)2  (2.38) 

may  be  calculated  and  used  to  estimate  Var[Ys]: 

Var[Yg]  =*  I  W2S*  /nk  .  (2.39) 

k=*l  K  K 

If  Yk  is  normally  distributed  with  mean  yyk  and  variance 
2 

°Yk'  then  an  approximate  100(1  -  6)%  confidence  interval 
for  Uy  j-3  given  by 

*3  1  t1_<5/2( ne  d.f.)  '  (Var[Ysl)1/2  (2.40) 

where  ne  is  given  by  [WELC56] 

L 

ne  ■  ($ar[Ya])2/  2  (wkSyk)4/In2(n  -  1)3  . 

k-1  k  k  (2.41) 

Stratification  is  generally  difficult  to  apply 


in  simulation.  While  many  methods  for  large  variance 
reductions  have  been  devised  for  specific  cases,  there 
are  no  general  techniques  for  queueing  simulations. 
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2.2.2  Control  Variates 

Control  variates  are  random  variables  with  known 
expectations  which  have  a  strong  linear  correlation  with 
some  response  variable  of  interest.  To  apply  control 
variables  to  a  simulation,  a  new  estimator  is  formed 
which  is  the  original  estimator  plus  a  linear 
combination  of  the  control  variables. 

Let  C  be  a  column  vector  of  Q  control  variables 

•St 

C  *  [Cl,  ...,  Cq]T  (2.42) 

with  expectation  vector  y^.  A  controlled  estimator  of 
the  simulation  response  Y  is  given  by 

Y(a)  *  Y  -  aT(c  -  yC)  .  (2.43) 

Y(a)  is  an  unbiased  estimator  of  yy  for  any  fixed  vector 
a  of  control  coefficients,  and  its  variance  is  minimized 
with  the  optimal  control  coefficient  vector 

®o  *  ^ C^2YC  (2.44) 

where  Zq  is  the  covariance  matrix  of  C  and  aye  is  the 
column  vector  of  covariances  between  Y  and  the  com¬ 
ponents  of  C.  The  minimum  variance  obtained  with  ao  is 
given  by 
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Var[Y(a0) ]  -  (1  -  R*c)Var[Y] 


(2.45) 


where 

R^c  *  ?YC?C1?YC/^Var^Y^  (2.46) 

is  the  square  of  the  coefficient  of  multiple  correlation 
between  Y  and  C. 

Under  ideal  conditions , 

Y(a0)  *  My  +  e  (2.47) 

where  e  is  an  irreducible  error  term  with  expected  value 
zero.  Combining  (2.43)  and  (2.47),  we  have  the  standard 
linear  regression  model 

Y  -  yY  +  a*(C  -  yc)  +  e.  (2.48) 


If  we  make  K  >_  Q  +  1  independent  replications  of 
the  simulation  and  we  let  Yfc  and  Cfc  be  the  observed 
values  from  the  kth  run,  then  we  will  obtain  K  iid 
observations  of  the  random  vector. 


The  {?*}  have  mean  vector 


and  covariance  matrix 


Let 


and 


‘  2 

aY  ?YC 
?YC  Ic 


X  * 


s  - 


Ken  -  uA)  ...  (cqi  -  Uq) 
•  •  • 

•  •  • 

•  •  • 

1(cik  "  ^l)  •••  ^CQK  “  WQ) 

W 

?o 


(2.51) 


(2.52) 


(2.53) 


Equation  (2.48)  may  therefore  be  expressed  as 

Y  ■  Xg  +  E  .  (2.54) 

The  least  squares  estimator  for  8  is 

8  -  (XtX)-1xty  .  (2.55) 

Lavenberg  [LAVE78]  was  able  to  show  that  the  straight¬ 
forward  estimator 


obtained  by  inserting  sample  covariances  into  (2.44) 
coincides  with  (2.55).  The  least  squares  estimator  of 
yy  is  therefore  given  by 

Y(a0)  -  Y  -  a£<C  -  uc)  .  (2.57) 

This  method  of  applying  control  variables  has 
some  problems.  Since  aQ  is  an  estimate,  minimum 

A 

variance  for  Y(a0)  will  not  be  achieved.  The  dependence 

A 

between  a0  and  C  generally  produces  a  biased  estimator 

A 

of  yy  Additionally,  dependence  among  the  Yfc(a0)  pre- 
vents  us  from  constructing  a  confidence  interval  using 
the  t  statistic.  There  are  two  methods  available  for 
handling  these  problems. 

If  Z  has  a  multivariate  normal  distribution, 
then  the  conditional  distribution  of  Y  given  C  •  c  is 
univariate  normal  with 

ElYj  C  -  c]  -  my  +  5o(S  ~  JJc*  (2.58) 

and 

VarlYj  C  -  cl  «  (1  -  rJc>Oy  *  (2.59) 
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We  nay  now  construct  a  100(1  -  S)l  confidence  interval 
for  uy! 

_  A  1/2 

Y(a0)  t  (2.60) 

where  sn  is  the  row  1,  column  1  entry  of  (XTX)-1  . 
Lavenberg  [LAVB78]  has  provided  an  expression  for  the 
loss  in  variance  reduction  caused  by  estimating  a0: 

Var[Y(a0)]/Var[Y(a0)]  »  (K  -  2)/{K  -  Q  -  2) 

(2.61) 

when  K  >  Q  +  2. 

A  second  technique  for  developing  confidence 
intervals  which  is  unconstrained  by  the  distribution  of 
Z  is  based  on  jackknifing.  Let  Y^tao)  be  the  estimator 
of  the  same  form  as  Y(ac)  but  with  the  single  vector 
Zk  omitted.  The  pseudo-values  are 

9*  -  KY(a0)  -  (K  -  l)Y-fc(a0)  ,  1  <  k  <  K  (2.62) 

with  the  point  estimate  of  Uy  given  by 

K 

W  -  (1/K)  l  0v  .  (2.63) 

k-1 

The  sample  variance  of  the  pseudo-values  is 


♦  '+! 
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% 


♦ 

* 


.  K 

a :  -  [1/{K  -  1)J  E  (0k  -  0>2  •  (2.64) 

0  k-1 

Thus  an  approximate  100(1  -  5)%  confidence  interval  for 
uy  is 

A 

0  ±  tl-fi/2;K-l  •  ae/Ki/2  .  (2.65) 

A 

Since  Yta^)  is  the  least  squares  estimate  of  ]iy,  this 
technique  is  equivalent  to  jackknifing  the  point 
estimator  of  the  first  method. 

2.2.3  Importance  Sampling 

Importance  sampling  completely  replaces  the 
original  sampling  process  with  another  one.  To  correct 
any  distortion  which  may  arise,  the  observations  are 
weighted  so  that  their  weighted  average  still  gives  an 
unbiased  estimate  of  the  mean  of  the  original  process. 
This  technique  is  somewhat  akin  to  stratified  sampling 
in  that  the  sampling  process  is  changed  and  the 
observations  are  weighted  differently. 

Suppose  the  objective  of  the  simulation  is  to 
evaluate  the  integral 

0  *•  fn  g(x)f(x)dx  (2.66) 


#*  *  **• 
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i 

« 


where  f(x)  is  a  density  function.  The  crude  Monte  Carlo 
procedure  would  be  to  randomly  sample  n  values  of  a 
variate  X  having  the  density  f(x)  and  take 

*  n 

9  »  (l/n)  2  g(Xji)  (2.67) 

i-1 

as  an  estimator  of  6 .  Let  h(y)  be  another  density 
function.  Then 


9*  f d  Kg(y)f (y) )/h(£)Jh(£)d£  .  (2.68) 

Let 

g*(y)  *  (g(y)f (£) )/h(jr)  .  (2.69) 

Then 

E[g*(y)l  *  /Dg*(^)h(y)dy  *  0  .  (2.70) 

Thus,  if  we  randomly  sample  n  values  of  a  variate  Y 

"m 

having  the  density  h(y),  our  estimator  of  0  is 

9h  -  (l/n)  2  g*<*i)  .  (2.71) 

i«l 

Minimum  variance  of  the  sampling  estimator  8h  is 
achieved  when  (KAHM53] 

&*<£>  »  I)  g(£)|  t{p)/U^  g(u)|  f(u)du)  .  (2.72) 


This  implies,  however,  that  we  already  know  the  value  of 
8.  If  h(y)  is  chosen  to  closely  mimic  g(y)f(y),  a  large 
variance  reduction  may  be  obtained.  If  h(y)  is  chosen 
poorly,  a  variance  increase  may  be  observed. 

Importance  sampling  was  developed  for  use  in 
Monte  Carlo  work  and  has  given  excellent  results  there 
[KAHN53,  KLEI74] .  However,  the  technique  is  more  dif¬ 
ficult  to  apply  to  discrete-event  simulation  and  the 
results  are  frequently  less  favorable. 

Moy  [MOY65]  proposed  a  "standard"  type  of  new 
density  function  to  be  used  in  importance  sampling  for 
simulation.  This  "standard"  function  does  not  require 
extensive  system  analysis  in  order  to  select  a  density. 
Suppose  X  is  the  response  variable  from  a  simulation  run 
which  uses  a  sequence  of  (pseudo)  random  numbers 
{?!,  . ..,  rm}  of  length  m.  Given  that  m  is  fixed,  we 
have: 

X  -  g(ri,  ..,  rm)  -  g  (R)  .  (2.73) 

The  expected  value  of  X  is  given  by 

E (Xl  *  /g  *»*  /q  g (*i,  ra)f (ri,  •••,  Em) 

•(dri  ...  dr*) 

«/  g(R)f (R)dR  . 

jta  -  •  - 


(2.74) 
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Since  the  random  numbers  are  independent*  we  have 


f ( r^ ,  • i < i  rm)  *  f 1(^1)  f u( rm)  ,  (2.75) 


where  f j  is  the  jth  marginal  density.  Th^se  marginals 
must  also  be  identical.  Let  fj(rj)  *  s(rj).  Then 


f <R)  *  fi(ri)  ...  fm(rra) 


»  s(ri)  . . .  s(rm)  .  (2.76) 


Since  we  are  dealing  with  random  numbers  uniformly 
distributed  over  [0,1],  we  have 


1,  0  <  r  j  <  1 

s(r j )  «  .  (2.77) 

0,  otherwise 


Therefore, 


II,  0  <  rj  <  1  for  all  j 

(2.78) 

0,  otherwise 


The  expectation  of  X  then  is 


EtXl  -  /  g(R)f (R)dR 

jia  -  ~ 

1  1 

*  /  ...  /  g(R)dR 

0  0  *  ~ 


1  1 

-  /  ...  /  (g(R)/h(R) ]h(R)dR  .  (2.79) 

0  0  *  *  *  ** 
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If  h(R)  is  a  joint  density  function,  we  may  sample  the 
input  vector  R  from  h(R)  and  use  the  importance  sampling 
estimator 

g*(R)  *  g (R)/h(R)  .  (2.80) 

Kleijnen  [KLE174]  uses  the  term  "importance  numbers"  to 
distinguish  between  stochastic  variables  with  an  alter¬ 
native  density  h(R)  and  random  numbers  having  the  stan¬ 
dard  density  (2.78).  The  variance  of  the  importance 
sampling  estimator  is 

Var lg*(R) ]  «/  m  (g2(R)f 2(r) )/h2(R) ]h(R)dR 

**  jIU  — 

-  (E [X] )2  .  (2.81) 

Selection  of  h(R)  is  again  critical,  and  Moy  recommended 
the  sampling  density 

|[in(  )/(o-l) )mexp{)ln(a) *Er j} 

0  £  r  j  £  1  £or  a11  3  (2.82) 

0  ,  otherwise 

where  a  is  a  parameter  which  he  estimated  by  using  a 
numerical  technique  to  solve  a  sample  version  of  the 
equation 
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■jj-  Var[g*(R;a)l  -  0  .  (2.83) 

dot  - 

In  M/M/1  queueing  systems,  Moy  obtained  variance 
reductions  of  33%  to  54%  using  this  approach.  He  also 
found  that  the  reduction  was  fairly  insensitive  to  the 
estimate  of  a.  He  also  showed  that  in  a  variety  of 

A 

systems  a  remained  fairly  constant  with  a  value  of  1.12. 
In  more  complex  systems,  Moy's  approach  produced 
variance  increases. 

2.2.4  Antithetic  Variates 

The  antithetic  variate  method  is  used  to  create 
negative  correlation  between  paired  observations  of  a 
response  variable.  One  observation  is  generated  from  a 
sequence  of  input  random  numbers  and  the  second  obser- 
vat ion  uses  the  corresponding  complimentary  sequence  of 
random  numbers. 

Suppose  y  is  the  mean  response  of  the  system. 

If  we  make  a  pair  of  simulation  runs  using  random  num¬ 
bers  {  rjj  on  the  first  run  and  {  (1  -  rj,)}  on  the  second 
in  order  to  obtain  the  responses  Xj,  and  X2,  then  we 
estimate  y  by 

» 

*  -  ^(Xx  4-  X2>  . 


(2.84) 


The  variance  of  X  is  given  by 

VartX]  -  i  Var [X^]  +  ^  Var  [X2]  +  j  Cov(X1,X2)  . 

(2.85) 

Thus  Var[X]  will  be  less  than  the  mean  of  two 
independent  replications  if  Cov(Xi,X2)  is  negative. 

Suppose  X  depends  on  Y  and  that  X  is  a 
monotonically  increasing  function  g^  of  Y.  If  Y  is 
generated  from  the  random  number  r  by  Y  ■  P”1(r),  then  Y 
is  a  monotonically  increasing  function  g2  of  r.  This 
implies  that  X  is  a  monotonically  increasing  function 
93  ■  0  g2  of  r.  Thus  a  high  value  of  r  yields  a  high 

value  of  X,  but  its  antithetic  partner  gives  a  low  value 
of  X.  Therefore,  intuitively,  g3(r)  and  g3(l  -  r)  are 
negatively  correlated. 

In  simulating  a  single-channel  queueing  system, 
Pritsker  [PRIT79]  points  out  that  rather  than  using  the 
actual  antithetic  pairs,  the  same  effect  is  achieved  by 
switching  the  random  number  streams  used  by  the  inter¬ 
arrival  and  service  time  processes.  This  results  from 
the  fact  that  long  service  times  increase  traffic 
intensity  while  long  interarrival  times  will  decrease 
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Antithetic  variates  should  be  used  with  some 
caution.  If  the  response  variable  X  has  a  symmetric 
distribution  and  the  inverse  transformation  method  is 
used  to  sample  X,  then  the  correlation  between  anti¬ 
thetic  partners  is  -1  and  a  100%  variance  reduction  is 
realized.  For  non-symmetric  distributions,  the  negative 
correlation  is  somewhat  weaker,  giving  smaller  reduc¬ 
tions.  In  some  cases,  the  use  of  antithetic  variates 
can  produce  a  variance  increase. 

2.2.5  Common  Random  Numbers 

Common  random  numbers  are  employed  when  the 
simulator  is  studying  more  than  one  system  and  needs  to 
choose  among  them.  The  responses  of  interest  are  not 
the  individual  system  responses,  but  the  differences 
between  them.  In  order  to  determine  these  variations, 
we  attempt  to  run  the  systems  under  the  same  conditions. 
To  accomplish  this,  the  same  initial  conditions  should 
be  used,  and  the  same  random  number  streams  should  be 
used  to  drive  similar  input  processes  (service  times, 
interarrival  times,  etc.).  This  usage  of  common  random 
numbers  results  in  positive  correlation  between  alter¬ 
native  system  responses,  say  Xi  and  X2.  The  variance  of 
their  difference  is  then  given  by 
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VarlXl  -  X2]  -  Var[Xi]  +  Var[X2] 

-  2Cov(Xi,X2)  .  (2.86) 

Compared  to  the  standard  approach  of  using  two 
independent  runs,  a  variance  reduction  will  be  observed 
if  a  positive  correlation  exists  between  X^  and  X2.  If 
large  values  of  r  result  in  large  values  of  Xj.  and  X2, 
then  a  variance  reduction  will  occur. 

In  these  simulations,  the  technique  requires 
that  we  operate  under  similar  conditions.  In  addition 
to  the  same  initial  conditions,  this  implies  that  the 
same  random  number  should  be  used  in  each  simulation  at 
the  same  juncture  in  the  operation  of  each  process. 
Kleijnen  [KLEI74]  states  that  this  "synchronization" 
could  be  maintained  more  easily  if  each  stochastic  input 
variable  has  its  own  random  number  stream.  Garman 
[GARM71]  suggests  using  a  single  random  number  stream 
but  discarding  some  numbers  to  maintain  synchronization. 

Large  variance  reductions  may  be  obtained  using 
this  method.  (Kleijnen  [KLEI74J  cites  reductions  as 
great  as  90  percent).  This  technique  is  advantageous  in 
that  additional  programming  and  run  time  are  minimal. 
This  method  has  been  used  in  conjunction  with  antithetic 
variates  with  mixed  results. 
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2.2.6  Application  of  Conservation  Equations 

For  regenerative  queueing  simulations,  the  use 
of  conservation  equations  to  estimate  steady-state  para¬ 
meters  can  produce  significant  variance  reductions. 
Carson  and  Law  [CARS80]  demonstrated  their  use  for 
Gl/G/s  queues.  In  the  GI/G/s  simulation,  values  of 
interest  include  the  mean  delay  in  queue  y^,  the  average 
number  in  queue  y^,  the  mean  sojourn  time  in  system  y^, 
the  average  number  in  system  uq,  and  the  average  amount 
of  work  in  system,  y^*.  Let  X  be  the  arrival  rate  and  y 
be  the  service  rate.  We  also  define: 

Y(l)  »  total  waiting  time  of  customers  in  the 
ith  cycle, 

<*i  ■  number  of  customers  served  in  the  ith 
cycle, 

£i  »  length  of  the  ith  cycle, 
ni  ■  total  service  time  of  customers  in  the 
ith  cycle, 

*J2>  .  yjl)  ♦  v 

and  Y|3)  ■  total  work  in  the  ith  cycle. 

Let  7(3.),  a,  %,  rw  7(2),  and  7(3)  be  their  respective 
sample  means  computed  over  n  cycles.  For  direct 
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simulation,  we  use  the  following  estimators  as  discussed 


in  section  2.1s 

0*  -  y<1)/«  ,  (2.87) 
UL  -  Y<U/T  ,  (2.88) 
VD  ■  Y(2>/o  ,  (2.89) 
UQ  -  Y<2 )/Z  ,  (2.90) 
GW*  -  Y<3)/C  .  (2.91) 


The  parameters  uw,  Ud#  Vq»  and  Uw*  are  related  by 
the  conservation  equations  for  the  system: 


UL  -  luw  ,  (2.92) 

Uo  ■  WW  +  (1/V)  »  (2.93) 

UQ  •  ^Ud  •  (2.94) 

and  UW*  -  (X/V)UW  ♦  (X/2)E[v2]  ,  (2.95) 


where  v  is  the  service  time  variate.  For  any  simulation 
of  a  GI/G/a  queue,  X,  y,  and  E(v2]  are  known.  Thus,  any 
estimate  of  any  of  the  five  parameters  may  be  used  to 
estimate  any  of  the  others.  Carson  and  Law  (CARS80] 
showed  that  the  most  efficient  estimator  of  uw  i*  w*# 
the  direct  simulation  estimator;  for  other  parameters. 


the  conservation  equations  were  shown  to  give  more 
efficient  estimators* 


ML  »  *MW  , 

(2.96) 

Md  *  My  +  (1/y)  , 

(2.97) 

Mq  *  Xyw  +  (X/y)  , 

(2.98) 

and 

yw*  -  ( X/y )yw  +•  (X/2)E[v2]  . 

(2.99) 

Using  this  scheme  with  1,  2,  and  4  servers, 
Carson  and  Law  [CARS80]  obtained  0  to  99%  variance 
reductions.  An  additional  advantage  of  this  method  is 
the  computational  savings.  Only  yj1*  and  cu  need  to  be 

a 

collected  and  used  to  compute 

2.3  Use  of  Control  Variables  in  Regenerative  Analysis 

Control  variables  have  been  applied  to  a  variety 
of  queueing  networks  with  varying  degrees  of  success. 

For  some  systems,  many  effective  control  variables  may 
be  found  and  easily  implemented.  In  other  cases, 
however,  it  may  prove  extremely  difficult  to  find  even 
one  such  variable  and  to  implement  it  if  one  is  located. 

Lavenberg  at  al.  [LAVE78]  explored  some 
controls  for  closed  queueing  networks.  These  networks 
have  a  finite  number  S  of  interconnected  service  centers 
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each  having  a  single  or  multiple  server  queue.  There 
are  a  finite  number  of  customers  N  which  circulate  among 
the  centers.  While  multiple  customer  types  were  con¬ 
sidered,  for  simplicity  w^  will  only  examine  the  single¬ 
type  case.  This  still  leaves  a  broad  class  of  networks 
since  there  are  no  restrictions  placed  upon  the  service 
time  distribution,  the  queueing  discipline,  or  the  queue 
capacity  at  each  station.  Some  assumptions  are  made: 

1.  The  sequence  of  entries  to  the  service  centers 
forms  an  irreducible  Markov  chain  with  a  state 
space  contained  in  the  set  of  service  centers 
and  a  fixed  transition  probability  matrix 

?  •  IPijl. 

2.  The  sequence  of  service  times  at  each  center  i 
is  a  sequence  of  iid  non-negative  random 
variables  distributed  as  the  random  variable  Tj.. 

3.  The  above  sequences  are  mutually  independent. 
From  pij's  the  long-run  relative  frequency 

with  which  a  customer  visits  center  i  may  be 
calculated  in  the  usual  way: 

VP  ■  ir  and  &i  •  1  •  (2.100) 

Lavenberg  et  al.  defined  an  event  to  be  the 


departure  of  a  customer  from  a  service  center,  end  the 
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customer's  service  time  is  associated  with  that  event. 
Let 


ejjM)  ■  number  of  events  associated  with  service 
center  i  in  the  first  M  events. 


and 

Wi(M)  *  sum  of  the  service  times  for  the  events 
included  above. 

Using  these,  Lavenberg  et  al.  developed  the  following 
control  variables: 

1.  W£(M)/e|{M),  a  service  time  variable  which 
represents  the  sample  average  service  time 
observed  at  center  i; 

2.  ejjMJ/M,  a  flow  variable  which  represents  the 
portion  of  events  occurring  at  center  i;  and 

3.  Wi(M)/M,  a  work  variable  which  represents  the 
ratio  of  completed  work  at  station  i  to  the 
number  of  events. 

Lavenberg  et  al.  proved  that 

lim  E(Wi(M)/ei(M>]  -  BCTiJ  ,  (2.101) 

M-*» 

lim  E[ei(N)/N]  -  *i  ,  (2.102) 

M-*» 


(2.103) 


lim  E[Wi(M)/M]  -  IfiElTil  . 

The  following  response  characteristics  were 
examined: 

wts  -  average  waiting  time  at  station  s, 

1  <  a  <  S  , 

A  “  steady-state  service  completion  rate,  and 
rt  *  average  time  for  a  customer  departing  a 
given  station  to  return  to  that  station. 

After  applying  the  control  variables  to  estimate  the 
above  performance  measures,  Lavenberg  et  al.  found  that 
the  work  variables  produced  larger  variance  reductions 
than  the  service  time  or  flow  variables  individually  or 
together.  Using  all  of  the  work  variables  defined  on 
the  network,  they  were  able  to  obtain  confidence  inter¬ 
val  coverage  comparable  to  the  uncontrolled  estimates 
based  on  the  Student-t  statistic.  Comparing  the  results 
obtained  from  multinormal  regression  theory  with  the 
results  based  on  jackknifing  regression,  they  found  that 
the  intervals  produced  by  the  first  method  were  signifi¬ 
cantly  narrower  and  required  much  less  computation  to 
obtain.  They  further  point  out  that  work  variables  give 


a  substantial  variance  reduction  provided  that  the  loss 
factor  (K  -  2)/(K  -  Q  -  2)  discussed  in  S  2.2.2  is  not 
too  large. 

Lavenberg  et  al.  [LAVE79]  studied  multiple 
control  variables  applied  to  regenerative  simulation, 
with  emphasis  on  obtaining  valid  confidence  intervals. 
The  GI/G/1  queue  and  several  central-server  models  were 
studied.  Under  uncontrolled  regenerative  simulation,  a 
steady-state  system  parameter  r  is  expressed  as  a  ratio 
of  expected  values 

r  -  Em/Elal  ,  (2.104) 

where  Y  and  a  are  appropriate  random  variables 
accumulated  over  a  single  regenerative  cycle.  (This  is 
the  standard  regenerative  notation  established  in  sec¬ 
tion  2.1.)  The  system  is  simulated  for  a  prespecified 
number  of  cycles  or  for  a  given  amount  of  simulated 
time,  and  the  pairs  {(*i,ai)f  are  collected.  Frequently 
several  long-run  average  system  parameters  cO*), 

1  <  k  <  R,  have  Known  values  and  can  be  expressed  as 

c(JO  -  EIY<*>]/Ela<*>j  ,  (2.105) 

where  aUO  end  y(*)  are  auxiliary  variables  defined  with 
respect  to  a  single  tour.  Let 


45 


« 


« 

* 


-s 

< 


_  n 

o  -  (1/n)  Z  oi i  ,  (2.106) 

i-1 

n 

Y  -  U/n)  Z  Y*  •  (2.107) 

i-1 

r  -  Y/a  ,  (2.108) 

Z|k)  «  yjk)  -  c(Jt)cinc>  ,  (2.109) 

n 

T(fc)  -  (1/n)  Z  Z(k)  ,  (2.110) 

i-1  1 

K 

and  r(a)  *  r  +  Z  avzW/a  ,  (2.111) 

k-1 


where  a  -  (a^,  a r).  Thus  r  is  the  standard 

regenerative  ratio  estimator  of  r,  r(a)  is  the  top* 

«v 

controlled  ratio  estimator,  and  (Z  (*)/»}  are  the  control 
variables.  Lavenberg  et  al.  established  that  r(a)  is 
strongly  consistent.  Let 


X 

o2(a)  -  VarlY  -  ra  ♦  Z  a^z^M  #  (2.112) 

k-1 


and 
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S2(a)  -  tl/(n  -  1)]  Z  {Yi  -  rai 

i-1 


+  Z  aic(Zi(J«>  -  ZOO)]2. 
k-1 


(2.113) 


Lavenberg  et  al.  state  that 


lira  S2(a)  *  o2(a)  with  probability  1  (2.114) 


and  that 


lira  Pr{/n  a(r(a) 

n-x» 


-  r)/S(a)  <  1}  »  *(z)  .  (2.115) 


From  this,  a  confidence  interval  may  be  constructed 


for  r. 


For  the  GI/G/1  queue,  the  response  variable 


studied  was  the  steady-state  mean  time  in  system,  li^. 


a  -  a(l)  ■  a(2)  ■  number  of  customers  served  in 
a  tour, 

Y  *  total  time  in  system  for  all  customers 
served  in  a  tour, 
y(l)  *  duration  of  a  tour,  and 
y(2)  *  duration  of  a  busy  period. 


If  the  arrival  rate  is  X  and  the  service  rate  is  y,  then 
E[Y(1.)]  -  E  [a]/A  (2.116) 

and 

E[Y(2)J  «  E[ct]/y  .  (2.117) 

The  control  variables  selected  were 

T(k)/a  -  [?(k)  -  C(k)a(k)]/S'(k) 

-  (Y<*>/a)  -  c(lc),  k  -  1,  2  ,  (2.118) 

where  c(D  »  1/X  and  c(2)  »  1/y. 

Lavenberg  et  al.  used  the  controls  with  the 
estimated  optimal  control  coefficient  a0  determined  by 
two  methods.  In  the  first,  or  dependent,  method  aQ  is 
estimated  from  the  first  m  of  the  n  (m  <  n)  tours  used 
to  construct  the  confidence  interval  for  r.  In  the 
second,  or  independent,  method  aQ  is  estimated  from  m 
tours  that  are  statistically  independent  of  the  n  tours. 
Independent  estimation  of  aQ  produced  smaller  variance 
reductions  than  were  obtained  with  dependent  estimation 
of  Sq.  Confidence  interval  coverage  was  maintained  at 
nominal  levels  with  independent  estimation  but  was 
significantly  less  with  dependent  estimation. 
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Xglehart  and  Lewis  [IGLE79]  also  proposed  a 
series  of  control  variables  for  regenerative  simulation. 
They  concentrated  their  efforts  on  the  estimation  of  the 
steady-state  mean  waiting  time  uw  in  the  GI/G/1  queue. 
Iglehart  and  Lewis  looked  for  controls  which  are  highly 
correlated  with  Yx  -  rctx  and  have  the  form  Cl  »  Di  - 
al/u,  where  U  is  the  service  rate.  Di  is  an  attempt  to 
mimic  Yj.  Using  the  recursive  relationship  (21) 
relating  successive  waiting  times,  they  suggested  the 
following  alternatives  for  Di* 


Da> 


*0  “  0  ,  -  1 
Wq  ♦  WX  ,  ox  >  2  ? 


D<2) 

1 


»0  *  0  ,  81  •  1 
W0  +  Wi  ,  -  2 

Wq  +  Wi  +  Xj  f«l  >.  3 


D<3)  >  W2  ; 


and 


D(4) 

1 


Wo  »  0  ,  <*1  «  1 
Wo  ♦  Wx  >«i  »  2 
Wq  ♦  Wi  +■  W2  $  ®1  >,  3 


(2.119) 


(2.120) 


(2.121) 


(2.122) 


*• 
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Thus,  the  controls  are  Ci d)  ■  Did)  -  ai/y  ,  i  * 

1,  2,  3,  4.  Iglehart  and  Lewis  point  out  that  it  is 

* 

0  more  difficult  to  calculate  E[Did)]  as  i  increases,  but 

that  Did)  more  closely  mimics  Yi  as  i  increases. 

In  testing  their  controls  for  an  M/M/1  queue, 
they  found  that  Ci<2),  and  Ci<3),  and  Ci<4)  performed 
much  better  than  Cid);  but  that  Ci<3)  and  Cid)  gave 
little  improvement  over  the  much  simpler  Ci<2).  There¬ 
fore,  they  selected  Ci<2)  as  the  most  desirable  control. 
Using  Ci<2)r  they  were  able  to  obtain  an  average  50% 
variance  reduction  in  the  estimation  of  mean  waiting 
time  over  a  wide  range  of  traffic  intensities. 

Wilson  [WILS79]  also  investigated  potential 
control  variables  for  regenerative  simulation.  The  net¬ 
works  he  considered  were  assumed  to  possess  finite 
expected  cycle  lengths  and  finite  expected  customer 
counts  at  each  station  for  each  cycle.  The  input  pro¬ 
cess  for  each  station  i  {X^d)  :  k  >  1)  (i.e.,  observed 
interarrival  times  or  service  times)  is  assumed  to  have 
a  known  distribution  wish  mean  and  variance  o*. 

Let 

gd,t)  *  number  of  service  times  that  are 
sampled  at  station  i  in  the  time 
period  [0,t), 


** 


so 


n(i,t)  -  number  of  service  times  that  are 

completed  at  station  i  in  the  time 
period  [0,t], 


n  (t)  -  Z  n( i,t) . 
i 


For  each  station  i,  Wilson  proposed  three  "standardized" 
control  variables.  The  first  of  these  is  given  by 

. . .  —1/2^  ^ 

Ci(U(t)  -  [g( i,t) ]  z  (Xk(i)  -  Ui)/oi  . 

k-1  (2.123) 


Wilson  was  able  to  show  that 


Ci(1>(t)  — ►  N( 0 , 1 ) 


(2.124) 


at  each  station  i.  His  second  control  is  a  standardized 
work  variable  for  Lavenberg's  closed  networks: 


Ci<2)(t)  ■  {  tn(irt)  J1/2/[n*(t)iri]} 


n(i,t) 

Z  (Xk(i)  -  yi)/ci  , 
k-l 


(2.125) 


where  w*  is  the  long-run  relative  frequency  with  which  a 
customer  visits  station  i.  For  such  a  closed  network. 


*  * 
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Ci ( 2 ) ( t )  ►  N (0 , 1 )  (2.126) 

t-H» 

* 

« 

at  each  station  i.  C^2)  takes  into  account  information 
about  network  traffic  flow  as  well  as  the  input  pro¬ 
cesses.  For  networks  in  which  some  subset  Kj  of  the 
stations  operate  independently  of  one  another ,  Wilson 
suggested  a  weighted  sum  of  control  variables: 

Cj  (3 )  ( t)  =»  E  wjiCi^Mt)  ,  (2.127) 

ieKj 

where  the  weights  wj^  are  arbitrarily  selected  by  the 
simulator.  This  control  variable  may  be  useful  in 
reducing  the  total  number  of  control  coefficients  which 
need  to  be  estimated  for  large  networks.  It  can  be 
shown  that 

o£ 

Cj(3>(t)  — ►  N{0,  E  wji2)  .  (2.128) 

fc*«  ieKj 

The  asymptotic  stability  of  these  standardized  work 
variables  allow  the  simulator  to  apply  either  poststra- 
tified  sampling  techniques  or  control  variate  analysis 
in  conjunction  with  regenerative  or  replication 
analysis. 

Wilson  was  able  to  obtain  variance  reductions 
ranging  from  30%  to  90%  in  a  variety  of  networks  when 


these  standardized  control  variables  were  applied  to  the 
numerator  of  classical  regenerative  ratio  estimators. 

At  the  same  time,  however,  he  encountered  significant 
estimator  bias  and  underestimation  of  the  variance 
resulting  in  loss  of  confidence  interval  coverage. 

These  difficulties  are  present  in  controlled  regenera¬ 
tive  analysis  regardless  of  the  controls  selected. 

Wilson  found  that  batching  the  cycles  resulted  in 
improved  coverage  while  still  giving  large  variance 
reductions.  In  addition,  he  provided  a  procedure  to 
determine  a  sampling  size  large  enough  to  guarantee 
adequate  coverage. 
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PROCEDURES  FOR  IMPROVING  BIAS  AND  CONFIDENCE  INTERVAL 
COVERAGE  PROBLEMS  IN  CONTROLLED  REGENERATIVE  ANALYSIS 

This  chapter  presents  an  analysis  of  the 
problems  encountered  using  control  variables  in  the 
numerator  of  regenerative- type  ratio  estimators.  Two 
types  of  control  variables  appear  in  the  literature 
[LAVE79] :  (i)  a  variable  which  is  an  estimator  of  a 

known  quantity  defined  with  respect  to  a  second  stochas¬ 
tic  system  where  the  known  quantity  is  an  approximation 
to  the  response  variable  of  interest;  Cii)  a  variable 
that  is  an  estimator  of  a  known  quantity  which 
defined  by  known  parameters  of  the  system  being  J. -ectly 
simulated.  We  have  restricted  our  research  to  the 
second#  or  concomitant#  control  variables.  A  class  of 
controls  is  selected  and  a  technique  is  developed  that 
improves  the  performance  of  those  control  variables. 

In  addition#  this  chapter  includes  a  theoretical 
development  for  a  multivariate  test  for  normality. 

Tests  which  have  been  proposed  in  other  research  are 
discussed  to  indicate  the  motivation  for  developing 
another  technique. 
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3.1  Concomitant  Control  Variates 

3.1.1  Selection  of  Controls 

A  wide  variety  of  control  variables  have  been 
proposed  in  past  research.  In  choosing  a  class  of 
control  variables,  it  is  necessary  to  establish  some 
criteria  for  comparison.  In  using  any  variance  reduc¬ 
tion  technique,  ease  of  implementation  is  highly  impor¬ 
tant.  Methods  requiring  less  computation  for  the 
machine  and/or  the  simulator  are  preferable.  It  is  also 
desirable  to  find  controls  which  are  effective  for  a 
variety  of  experimental  systems.  For  these  reasons, 
only  concomitant  (or  internal)  control  variables  were 
considered  for  use  in  this  research.  Lavenberg  et  al. 
[LAVE78]  offered  an  effective  set  of  controls  for  closed 
queueing  systems  with  the  work  variables  they  proposed. 
While  these  controls  are  easily  ing>lemented,  they  are 
asymptotically  unstable  and  they  cannot  be  extended  to 
open  or  mixed  systems  [WILS79] . 

Lavenberg  et  al.  [LAVE79]  examined  multiple 
control  variables  applied  to  the  numerator  of  a  regener¬ 
ative  ratio  estimator.  These  controls  were  applicable 
to  both  open  and  closed  systems.  These  variables 
yielded  substantial  variance  reductions  in  the  models  to 
which  they  were  applied.  However,  their  controls 


produced  a  significant  degradation  in  confidence 
interval  coverage. 

Iglehart  and  Lewis  [IGLE79)  proposed  another 
class  of  controls  which  uses  a  recursive  relationship 
for  successive  waiting  times  that  is  only  applicable  to 
the  GI/G/s  queue.  While  their  controls  are  effective, 
derivation  of  the  expected  values  of  these  variables 
presents  a  formidable  task. 

Wilson  [WILS79]  offered  yet  another  type  of 
concomitant  control  variable.  The  only  restrictions  he 
levied  on  the  type  of  systems  to  be  considered  were  that 
they  possess  finite  expected  cycle  length  and  finite 
expected  visit  counts  per  cycle  at  each  station.  His 
controls  performed  effectively,  but  he  did  encounter 
some  loss  of  confidence  interval  coverage.  This  loss 
was  not,  in  general,  as  large  as  that  experienced  by 
Lavenberg  et  al.  [LAVE79]  under  their  complete  dependent 
estimation  scheme. 

With  this  collection  of  choices  of  concomitant 
control  variables,  we  selected  Wilson's  "standardised 
service-time  variates"  to  use  in  our  research.  While 
these  controls  produced  large  efficiency  increases  when 
applied  to  the  numerator  of  a  regenerative  ratio  estima¬ 
tor,  they  did  not  perform  as  well  when  applied  to  the 
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denominator  of  such  ratios.  For  this  reason,  we  propose 
the  use  of  a  "standardized  flow  variate"  to  control  the 
denominator.  For  example,  if  customers  move  from  sta¬ 
tion  1  to  station  2  with  probability  Pi2*  or  elsewhere 
with  probability  qi2  ■  1  -  P12,  then  each  branching 
operation  may  be  regarded  as  a  Bernoulli  trial  with  mean 
P12  and  variance  Pi2<3l2  (here  station  2  is  regarded  as 
a  "success” . ) 


Define 

h(j,t)  =  number  of  customers  which  arrive  at 

branch  j  in  the  simulated  time  period 
[0,t]  (3.1) 


Ifc(j)  -  the  indicator  function  for  a  "success" 
at  branch  j  for  customer  k  (3.2) 


Pj  =  probability  of  a  "success"  at 
branch  j 


(3.3) 


s  1  -  Pj 


(3.4) 


Thus,  Xk(j)  has  mean  Wj  -  pj  and  variance  o*  *  pjqj. 

J 

The  "standardized  flow  variate"  is  defined  to  be 


Dj(t) 


ih( j ,t) ) 


b{Iv(j)  -  Wj)  /  Oj  (3.5) 
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This  variable  is  o£  the  same  form  as  Wilson's 
"standardized  service-time  variate" 


v.g(j#t) 

Cs(t)  -  {g(j,t)J  ***&  E  {Yk(j)  “  • 

3  3c-l  J  (3. 


6) 


which  was  defined  in  Chapter  2. 


3.1.2  Theoretical  Development  of  the  Joint  Distribution 
of  Standaradized  Variables 

Suppose  we  have  a  regenerative  system  with  QC 
stations  and  QD  branching  points.  Thus,  we  may 
construct  QC  controls  of  the  form  (3.6)  and  QD  controls 
of  the  form  (3.5). 

Let 


Q  -  QC  +  QD 


(3.7) 


a( j,t) 


g  ( j  r  t )  ,  l  <  j  <  qc 

h(j  -  QC,t)  ,QC+l<j<QC+QD 

(3.8) 


UjcCj)  - 


**(j)  *  l  <  j  <  QC 

I)c<j  -  QC)  ,  QC  +  1  <  j  <  QC  +  QD 

(3.9) 


Thus  the  jth  standardized  control  may  be  expressed  as 


,a( j,t) 

Aj(t)  -  {a( j, t)J  ^  {ok(j)  - 


10) 
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and  the  Q-dimensional  control  vector  is 

A(t)  -  [Ci<t),  ...»  CqcU),  Di(t),  ...»  DQD(t)]T 

-  (Ai(t),  ...»  Aq( t) )T  .  (3.11) 

We  wish  to  show  that  A(t)  converges  in  distribution  to  a 
multinomial  distribution 

A(t)  -  Nq(0,  E)  (3.12) 

with  null  mean  vector  and  a  correlation-type  covariance 
matrix  E  .  Relation  (3.12)  will  provide  a  theoretical 
foundation  for  applying  multinomial  regression  theory  to 
the  study  of  controlled  ratio  estimators. 

We  begin  with  the  observation  that  [RA073] 

Z  -  Nq(  y  2 /  5  z ) 

<«>  bTz~N(bTyz,  b*  E-b)  for  all  be  RQ  .  (3.13) 

We  have  assumed  that  our  regenerative  queueing  system 
has  a  finite  non-zero  asymptotic  sampling  rate 

a  -t  =  lim  a(j,t)/t  a.s.  (3.14) 

for  the  process  associated  with  control  j#  1  <  j  <  Q. 
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Define  the  partial  sums 


n 

Sn(j)  2  E  {Uk<j>  -  Uj}  /  o*  ,  l  <  j  <  q  (3.15) 
k«l  J  -  - 


and  let 

[  a-st]  =  greatest  integer  in  a^tf  1  <  j  *<  Q  . 

“(3.16) 

(The  use  of  brackets  to  denote  the  greatest  integer 
function  is  limited  to  this  section  only.)  Note  that 
for  any  t,  the  variates  tS{0^tj(j)  :  1  <  j  <  q)  are 
mutually  independent  and 

Zj(t)  =  S(  ajti(j)/[  ojt]^^  N(0,  1)  (3.17) 

by  the  central  limit  theorem.  Let 

Z(t)  -  (Zi(t),  ...»  ZQ(t)]T  .  (3.18) 

Define  the  characteristic  functions 

♦fct(u)  =  Etexpt  iuZjc(t)}  ]  (3.19) 

(Note  that  in  this  section  only,  we  reserve  the  symbol  i 
to  denote  /  -1. )  For  a  fixed  b  in  rQ  we  have 


$t(u)  »  B (exp  {iubTZ(t)>  ) 

Q 

■  it  (3.20) 

k-1 

since  the  {Zk(t)}  variates  are  independent.  In  view 
of  (3.17) ,  the  continuity  theorem  for  characteristic 
functions  [NEUT73]  implies  that 

lim  ♦^(u)  *  exp  ^”^2/2)  for  all  u  .  (3.21) 

t-**» 


Combining  (3.20)  and  (3.21),  we  have 


w  2 

lim  $t(u)  *  exp  {(-u2/?)  I  bfc  > 
t-M  k»l 


(3.22) 


which  is  the  characteristic  function  for  a  normal 

Q  2 

variate  with  mean  zero  and  variance  I  b.  .  Again 

k-1  K 

invoking  the  continuity  theorem,  we  have 


b?Z(t)  ♦  bTZ  for  all  be  rQ?  Z  -  Nn(0,  Iq)  (3.23) 


where  Iq  is  the  Q  x  Q  identity  matrix.  To  prove  (3.12), 
consider  the  dissection  formula  for  A^(t) 


sa<j,t)  _  St  ajtl  ([  ajtJ  ^  ^  S(  ojtl 

(a(  j  ft)  )  ^  [  Ojtl*.  \a(j,t)/  t  “jt]% 


t  Qjt}^ 


i  gjtl  i 
a( j,t) | 


(3.24) 


Note  that 


lim  {[  autl/a(  j,t)  }  -  1  ,  1  <  j  <  Q  .  (3.25) 

t-«» 


In  view  o£  (3.17)  we  have 


Wj(t) 


U(j,t)  J 


sl  ajt)  p 

1  *  - J-rr.-  0  (3.26) 

[  Ojt]^ 


by  Slutsky's  theorem  [BICK77].  Wilson[WILS79]  proved 


that 


r3,t,  ,  {;■«>»  -s' 1 a*L}.(l2!L}  !, 

J  (  (  o^tiw  )  (  a( jr t))  <3.: 


If  we  define  the  random  vectors 


W(t)  -  lWX(t),  ...,  WQ(t)]* 
V(t)  -  lVi(t),  ...,  VQ(t)l* 


then  relations  (3.26)  and  (3.27)  together  with  Slutsky's 
theorem  imply 
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P  P 

b^W(t)  -*■  0  t  b^V(t)  -*■  0  for  all  b  e  rQ*  (3.28) 
Now  the  dissection  formula  (3.24)  may  be  expressed  as 

A(t)  -  Z(t)  +  W(t)  +  v(t)? 
and  combining  this  with  (3.23)  and  (3.28),  we  have 

b^A( t)  *«  b^z  for  all  beRQ,  z  -  Nq(0,  Iq).  (3.29) 
Proposition  2c. 4  (xi)  of  [RA073]  finally  implies 

A(t)  -  Nq(0,  Iq)  . 

3.2  An  Examination  of  Bias  and  Confidence  Interval 
Coverage 

Previous  research  in  the  area  of  controlled 
regenerative  analysis  has  shown  that  large  variance 
reductions  may  be  achieved  through  the  use  of  con¬ 
comitant  control  variables  applied  to  the  nisoerator  of  a 
ratio  estimator.  Serious  problems  have,  however,  been 
encountered  in  their  use.  A  significant  level  of 
relative  bias  seems  to  be  inherent  in  the  top-controlled 
estimator  [WILS79].  Additionally,  it  appears  that  the 
confidence  intervals  obtained  frequently  fail  to  provide 
proper  coverage  probabilities  (LAVB79,  WILS79]. 
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3.2.1  Background . 

We  assume  that  the  queueing  system  to  be  studied 
possesses  the  regenerative  property,  that  the  I ID  cycle 
lengths  {<**,  k  £  1}  satisfy 


0  <  E[  ak]  < 


(3.30) 


and  that  Xk(i),  the  number  of  customers  served  at  sta¬ 
tion  i  during  cycle  k,  satisfies 


0  <  E [Xk( i) ]  <  *  for  all  i. 


(3.31) 


These  are  relatively  mild  constraints  upon  the  allowable 
queueing  systems. 

In  top-controlled  regeneration  analysis, 
response  variables  Y  and  X  and  a  vector  of  standardized 
control  variables  C  with  a  known  expectation  of  yc  *  0 
are  collected  for  each  of  the  n  simulated  tours.  The 
regenerative  ratio  r  =  yy/  yx  is  then  estimated  by 


r(  6  )  -  (Y  -  6  TC)/X 


where  Y,  X,  and  C  are  the  sample  means  over  the  n 


(3.32) 


cycles.  An  approximation  to  the  variance  of  r(  3  )  for 
large  samples  is  [IGLE79] 

V(r(  §  )]  -  tl/(n  t^)}  .  V[Y  -  rX  -  6  TC] .  (3.33) 

x  ~ 


64 


[ 

l  * 


« 


I 

i 

r 


The  variance  of  r(B)  is  minimized  with  the  optimal 
control  coefficient  vector 

§o  -  l p”1  <?(*  -  rX,  C)  ,  (3.34) 

where  £c  is  the  covariance  matrix  of  C  and  a(Y  -  rX,  C) 
is  the  column  vector 

[Cov(Y  -.rX,  CiH 

:  .  (3.35) 

Cov(Y  -  rX,  Cq)J 

The  optimal  control  coefficient  @Q  is  estimated  by 

b  -  E-l  o(Y  -  rxX,  C)  (3.36) 

using  the  corresponding  sample  covariances  from  the  n 
observed  cycles.  The  sample  variance  estimator  for  r(b) 
is  given  by 

A  A  n  __ 

V[r(b) ]  -  {1/ tn* (n-1) »X2J }  •  E  [Yj-r.Xj-bTtCj-C) )  2 

j-1  J  1  ~  (3.37) 

where  r^  *  Y/X.  Now  r(b)  is  asymptotically  normal 
(LAVE79] ,  and  therefore  a  100  (1  -  a)%  confidence 
interval  for  r  is  given  by 

r(b)  ±Zl-o/2*  Wr(b)  1  Vi  .  (3.38) 
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3.2.2  Theoretical  Examination  of  Bias  in  the  Top- 
Controlled  Estimator 

In  order  to  examine  the  bias  present  in  r(b), 
the  top-controlled  estimate  of  r,  let  us  first  discuss 
bias  as  related  to  the  classical  ratio  estimator 

ri  -  Y/X  (3.39) 

taken  over  the  same  n  cycles.  Let 

ri  =  Etrxl  .  (3.40) 

The  covariance  between  ?x  and  X  is  given  by 

Cov(ri,X)  m  E[ri  •  X]  -  E[rxl  •  E[X] 

-  E ( (Y/X)  •  X]  -  ri  yx 
■  Vy  -  ri  wx  .  (3.41) 

This  gives  us  the  relation 

(1/  yx)  •  Cov(ri,  X)  »  Uy/  Ux  -  ri 

*  r  -  rx  (3.42) 

Thus,  the  bias  in  rx  is  given  by 

B(ri)  =  rx  -  r  •  -  Cov(ri,  X)/  wx  .  (3.43) 
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which  is  of  the  same  form  as  (3.41).  Through  a  similar 
argument  we  obtain 

|  B[r(b)  1  |  -  {  1  R2  I  •  Setr(b)]  •  SE[X}  >  /  Ux  r 

(3.51) 

where  R2  is  the  coefficient  of  correlation  for  r(b)  and 
X.  This  yields  the  bound 

|B[r(b)]  |  /SE[r(b) ]  <  CV(X)  .  (3.52) 

3.2.3  A  Technique  for  Controlling  Bias 

Research  in  the  area  of  concomitant  control 
variables  has  been  confined  to  the  application  of 
controls  to  the  numerator  of  a  ratio.  Equations  (3.47) 
and  (3.52)  reveal  that  the  relative  bias  present  in  the 
classical  and  top-controlled  regenerative  ratios  primar¬ 
ily  relates  not  to  the  numerator  but  rather  the  denomi¬ 
nator  of  the  ratio  of  interest.  We  have,  therefore, 
chosen  to  turn  our  attention  to  the  development  of  a 
strategy  pertinent  to  the  denominator  which  would 
improve  the  efficiency  of  a  queueing  simulation  in  the 
same  manner  as  the  top-controlled  method  without  the 
bias  penalty. 

A  logical  path  of  exploration  is  to  consider 
applying  controls  to  the  denominator.  Let 


where  7  indicates  the  sample  mean  of  a  vector  of  stan¬ 
dardized  control  variables  D  with  known  expectation 
UD  *  0.  Applying  the  same  argument  used  earlier  we 
have 

Btr(  6  )]  »  -  Cov[r(  6  )  ,x  -  6^5]/  Mx 

-  -{R3  •  SE[r(  <S  )J  •  Se [X  -  /  Ux 

(3.54) 

where  R3  is  the  coefficient  of  correlation  between  r( 6  ) 
and  X  -  <$TP.  Thus,  a  bound  on  the  relative  bias  is 
given  by 

|  B(f(  6  )]  |  /SE [r(  6  )] 

-  SE [X  -  6TD]/(n^  Ux)  .  (3.55) 

As  will  be  seen  in  Chapter  V,  the  relative  bias  of  a 
ratio  estimator  is  a  major  cause  of  coverage  degradation 
in  the  associated  confidence  interval  estimator.  (See 
also  (COCH77J ,  pp.  12-16.)  This  last  expression  implies 
that  applying  controls  to  the  bottom  of  a  ratio  estima¬ 
tor  should  reduce  bias  if  we  are  able  to  find  controls 
which  are  strongly  correlated  with  X.  The  optimal 
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control  coefficient  vector  that  minimizes  the  bound 
(3.55)  on  the  relative  bias  of  r( 6  )  is  given  by 

«  O  •  S"1  2(X,  D)  »  (3.56) 

D 

where  Z d  denotes  the  covariance  matrix  of  the  random 
vector  D  and  cr(X,  D)  denotes  the  column  vector  of 
covariances  between  X  and  each  component  of  D.  This 
optimal  control  vector  is  estimated  by 

d  -  £-1  a (X,  D)  (3.57) 

~  D 

where  the  entries  of  E&  and  a(X,  D)  are  the 
corresponding  sample  covariances  computed  over  all 
observed  tours.  In  addition,  care  must  be  taken  that 
application  of  those  controls  does  not  increase  the 
correlation  between  the  rat  o  estimator  and  its 
denominator. 

3.2.4  A  Two-Stage  Procedure  for  Controlling  Bias  and 
Variance 

In  practice,  application  of  controls  to  the 
denominator  produced  strong  results  in  terms  of  bias 
correction.  These  controls  did,  however,  simultaneously 
result  in  large  variance  increases.  This  led  us  to 
attempt  a  two-stage  procedure.  Since  application  of 
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controls  to  the  numerator  gives  large  variance  reduc¬ 
tions  ,  we  propose  that  the  simulator  should  apply 
appropriate  controls  first  in  the  denominator  to  reduce 
the  bias  of  the  regenerative  estimator;  then  he  should 
apply  a  different  set  of  controls  to  the  numerator  of 
the  bottom-controlled  ratio  to  obtain  a  variance 
reduction. 

The  new  ratio  estimator  offered  by  this  research 
is 

r(  8,Y  )  *  (7  -  $TC)/(x  -  yTo)  (3.58) 

where  C  and  D  are  disjoint  sets  of  standardized 
controls.  The  variance  of  this  estimator  is  given  by 

V[r<  BfY  )]  5  V[Y  -  BTC  -  r(X  -  YTD)]/(n  y2). 

~  '  '  ~  ~  ~  x 

(3.59) 

We  first  wish  to  control  the  denominator  to  obtain  a 
reduction  in  bias.  This  requires  a  minimization  of 
V(x  -  YTp).  The  optimal  control  coefficient  is  there¬ 
fore  estimated  by 

d  -  E"1  a(x,  D)  .  (3.60) 

—  -vQ  —  — 


Let 


We  may  now  express  our  ratio  estimator  as 

r(  6»  d)  -  (Y  -  8T£)/X*  (3.62) 

which  has  approximate  variance 

V[r(  8,  d)]  -  V[(Y  -  6TC)  -  rX*]/n  u2).  (3.63) 

~  ~  rx 

The  vector  80  which  minimizes  (3.63)  is  estimated  by 

b  »  J-l  o(Y  -  riX*,  C)  .  (3.64) 

The  resulting  sample  variance  is  given  by 

V(f(b,  d)]  -  {  l/[n  •  (n  -  1)  •  (X  -  dT5)2] } 
n 

•  I  (Yj  -  bT(Cj  -  C)  -  ri[Xj  -  dT(Dj  -  D)]>  2 
j»l  ~  ~  (3.65) 

We  form  an  approximate  100  (1  -  a  )%  confidence  interval 
for  r  by 

r(b,  d)  +  Zx.a/2  -  (v[r(b,  d)]>  %  .  (3.66) 

The  relative  bias  structure  of  the  two~stage 
controlled  ratio  estimator  reveals  much  information  as 
to  the  settings  in  which  the  procedure  will  be 


beneficial.  For  the  first  time,  an  explicit  expression 
for  bias  is  given  which  will  allow  the  simulator  to 
determine  after  a  short  preliminary  run  whether  or  not 
the  controls  he  selects  will  result  in  bias  reduction. 
The  bias  of  our  estimate  is  defined  to  be 

B[r(b,  d)]  =  E(r(b,  d)]  -  yy/  jjx 

-  -  Cov[r(b,  d),  X  -  dTp)/  Vx.  (3.67) 

Now, 

Cov[r(b,  d),  X  -  dTD] 

-  R4  •  SE [r(b,  d)l  •  SE[X*1  ,  (3.68) 

where  R4  is  the  coefficient  of  correlation  between 

r(b,  d)  and  7*.  Therefore,  the  measure  of  relative  bias 

is 

|  Blr(b,  d)]  )  /SE[r(b,  d)] 

-  |  R4  |  •  SE[X*l/(ntt  Vx)  .  (3.69) 

Using  Lavenberg's  [LAVE78]  formula  for  loss  of 
efficiency  due  to  the  estimation  of  the  control  coef¬ 
ficients  over  n  regenerative  cycles,  we  obtain 


VfX*l  -  V[X  -  dTD} 

«  a  2  •  U  -  r2)  .  [(n  -  2)/{n  -  QD  -  2)), 
x  5  (3*70) 

where  we  assume  D  to  be  a  QD  dimensional  vector  and 
Rs  is  the  coe££icient  of  multiple  correlation  between  X 
and  D.  We  therefore  obtain  the  explicit  expression  for 
relative  bias 

|  B[r(b,  d) I  |  /SE(r(b,  d)l  » 

iox  •  Jr4|  •  (1-R|)*i  •  Un-2)/(n-QD-2)]^}/(n^ux) 

(3.71) 

This  last  expression  gives  several  clues  for  selecting 
controls  for  the  denominator.  Due  to  the  presence  of 
the  loss  factor  (n  -  2)/(n  -  QD  -  2 ) ,  keeping  the  number 
of  bottom  controls  to  a  minimum  is  desirable.  In  addi¬ 
tion,  a  strong  correlation  between  X  and  D  is  prefer¬ 
able,  but  it  should  be  accompanied  by  a  weak  correlation 
between  r(b,  d)  and  X*. 

To  gain  insight  into  the  behavior  of  the 
variance  of  the  two-stage  estimator,  we  will  now  develop 
an  explicit  expression  for  that  quantity.  This 
expression  will  again  give  the  simulator  information 
about  the  value  of  hia  controls. 
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i 

Let 

* 

i 

Mow, 

Yi*  «  Yi  -  b^Ci,  1  <  i  <  n 

(3.72) 

• 

V[r(b,  d)]  -  [l/(n  yx2)]  •  V [ (Y 

»*»  ■ 

-  bTC) 

-  r(X  -  dTD) ] 

•>» 

(3.73) 

and,  in 

havex 

terms  of  Y*  =  Y  -  bTC  and  X*  s 

X  -  dTD,  we 

M* 

V[(Y  -  bTC)  -  r(X  -  d*D)) 

•»  **  •# 

•  vlY*  -  rx*l 

-  V[Y*]  +  r2V[X*]  -  2r 

Cov(Y* ,  X*) 

»  °yU  -  Rg)Lr  +  r2  °*(1  -  Rj)Lb 

i 

-  2r [R7  -  a  (1  -  R^)^ 

*  •  Ojj 

!  * 

•  (1  -  r|)V5  .  lb^] 

(3.74) 

! 

1 

!  < 

whara  R$  and  R7  respectively  danota  tha 

correlation 

! 

* 

I 

i 

coefficients  batwaan  Y  and  C,  and  botwaon  Y*  and  X*. 

Moreover ,  we  Hava  tha  lost  coafficianta  for  tha  top  and 

* 

botton 

of  r(b,  d)i 
<* 

*  • 


« 


LB  *  (n  -  2)/(n  -  QD  -  2)  . 

Equation  (3.74)  reveals  that  there  are  several 
opportunities  for  variance  reduction  in  the  two-stage 
estimator.  Strong  correlations  between  Y  and  C  and 
X  and  D  will  reduce  the  overall  variance.  In  addition, 
a  strong  positive  correlation  between  Y*  and  X*  will 
result  in  variance  reduction,  but  a  variance  increase 
will  be  realized  if  they  are  negatively  correlated. 

We  will  point  out  here  that  results  in  S  3.1  are 
the  first  quantification  of  the  bias  found  in  these 
controlled  ratio  estimators.  Implications  for  the 


practitioner  will  be  summarized  in  Chapter  V. 

3.3  Development  of  a  Test  for  Multivariate  Normality 


The  class  of  control  variates  selected  for  use 
in  this  research  requires  that  they  be  from  a  multi¬ 
variate  normal  distribution.  This  required  the  batching 
of  observations  over  a  number  of  cycles  to  induce  a 
central-limit  effect.  In  order  to  reduce  the  total 
number  of  simulated  cycles  required,  we  wished  to  use  a 
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normality  test  which  would  indicate  the  minimum  required 
batching  size. 

* 

3.3.1  Analysis  o f  Previously  Proposed  Tests 

Hensler  et  al.  [HENS77]  developed  a  test  for 
multivariate  normality  requiring  successive  orthogonal 
transformations  of  the  observations.  They  base  their 
test  on  the  following  theorem  and  corollary: 

Theorem.  Let  Xi,  X2,  . ..,  £n  be  independent  k- 
variate  random  variables  with  covariance  matrix 
£  .  Let  A  *  (aij),  i,  j  *  1,  . ..,  n,  be  an  n  x  n 
orthogonal  matrix  over  Rl  such  that  at  least  two 
elements  in  the  first  row  are  not  zero.  Denote 
*iT  *  (Xiif  ...»  Xi]{)  i  i  ■  1,  ...»  n»  X  *  (Xij)» 
i  «  1»  . . . »  n,  j  -  1,  ...»  k  and  let 

X  -  AX 

where  Yi,T  .  (yilf  ...»  Yi*),  i  •  1,  ...»  n. 

Then  Y2,  X3»  •••»  Xn  are  (n  -  1)  iid  k-vartate 
normal,  N(0,  I  )»  random  variables  if  and  only 
if  Xi,  ...»  Xn  are  such  that  Xj  is  k-variate 
normal,  NUijBtYj.],  § ) ,  j  *  1,  . . . ,  n. 

Coilary.  Let  X],,  ...»  Xn  be  a  random  sample 
from  a  k-variate  distribution  with  mean  vector 


*  % 
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W  and  covartance  matrix  £  .  Define  an  (n  -  1) 
by  n  matrix  B  ■  (bj,j)  by 

1  -  b,  i  »  j,  i  «  1,  2,  n  -  1 

"  b'  i  ^  j#  i/  j  -  X,  2,  ...r  n  -  1 

“(1  +  /n)b, i  ■  1/  •••/  n  -  l;  j  a  q 

where  b  ■  l/[n  +  /"I n] .  Complete  the  matrix  B  to 
an  n  by  n  orthogonal  matrix  C.  Define  n  random 
variables  Yi,  Yn  by 

l  *  cx 

where  Y  -  (YijUnd  X  •  (Xij),  i  «  1,  n, 

j  -  lr  ...»  k.  Then  YX,  ...»  Yn  -  i  are  iid 
N(0,  Z  )  if  and  only  if  Xlr  xn  are  iid 

N(M,  E  ). 

This  corollary  allows  them  to  test  for  N (0,  Z  )  rather 
than  testing  for  N(  y,  E  ). 

<**  mt 

To  test  for  k-variate  normality  of  a  random 
sample  of  size  n,  the  corollary  must  first  be  employed 
to  reduce  the  problem  to  one  of  testing  for  a  sample  of 
size  n  -  1,  YtT  -  (Yilr  Yi*),  i  -  lf  n  -  1, 

from  a  tl(0,  E  )  distribution.  Hensler  et  al.  then 


1 

i 

j 


propose  testing  the  last  component  X(k)T  *  (Xik,  ••• 
Xn-i,]c)  for  univariate  normality  with  mean  0  and 
variance  oj^2.  Next,  they  test  the  conditional  (k  -  1) 
variate  normality  of  the  first  k  -  1  components  given 
the  last.  They  form  am  orthogonal  matrix  A  whose  first 
row  is 

xlk  x2k  ...  x(n  -  l)k 

s  s  s 

where 

n-1 

s2  *  E  Xik2 
1*1 

Let 

W  -  AX  . 

By  their  theorem,  this  reduces  the  problem  to  testing  W 
for  (k  -  l)-variate  normality.  The  procedure  is 
repeated  testing  the  last  component  for  univariate  nor¬ 
mality  and  forming  a  new  orthogonal  transformation. 

This  test  presents  several  problems.  First,  it 
omits  a  great  deal  information,  univariate  tests  for 
normality  are  performed  on  k  independent  sample  sizes  of 
n-1,  n  -  2,  ...,  n  -  k,  so  that  [k(k  +  l)/2]  obser¬ 
vations  are  lost.  Second,  the  practicality  of  forming 


*  * 
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T 


* 


0 


several  orthogonal  transformations  is  questionable. 

Given  large  sample  sizes,  this  could  prove  extremely 
expensive.  Third,  the  choice  of  the  test  for  univariate 
normality  is  left  to  the  tester.  Hensler  et  al.  offer 
no  information  on  the  power  of  this  test  but  merely 
state  that  it  will  do  no  better  than  the  univariate  test 
selected. 

Cox  and  Small  [COX78]  proposed  goodness  of  fit 
tests  for  multivariate  normality  based  upon  tests  of 
linearity  of  regression.  They  proposed  both  coordinate- 
dependent  and  invariate  approaches.  Due  to  the  theoret¬ 
ical  power  and  value  of  invariant  methods  in  multi¬ 
variate  analysis,  we  only  present  that  particular  test. 
The  basic  idea  of  their  method  is  to  find  the  pair  of 
variables  which  are  linear  combinations  of  the  original 
v?  k  \oles  that  will  result  in  maximum  curvature  when  one 
is  regressed  on  the  other.  The  amount  of  that  curvature 
is  taken  as  the  test  statistic. 

Suppose  Y  *  (Yi,  ...,  YV)T  is  a  standardized 
v-dimensional  variate  with  mean  0  and  covariance  matrix 
E  .  For  the  higher  moments,  write  for  r,  s,  t,  u  »  1, 
...»  v , 


E[YrY8Ytl 


y(r,s,t),  E[YrYaYtYu] 


- 


y(r,s,t,u) 
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Two  linear  combinations  X  3  a^Y  and  W  3  bTY  are  formed 
with  aT  Za  3  bT  Zb  3  1  so  that  X  and  W  have  mean  zero 
and  variance  one. 

Let  y  3  Yxw  be  the  least  squares  regression 
coefficient  of  X  on  w2,  adjusting  for  the  linear 
regression  of  X  on  W.  Using  the  orthogonal i zed  form 

X  3  8  W  +  Y  [W2  -  WE[W3]  -  1]  +  e  , 

Cox  and  Small  obtained  the  coefficient 

E[X  •  W2]  -  E[W3]  •  E[X  •  W] 

Y  =  - - - -  — — — .  ■  * 

XW  E[W4]  -  1  -  (E[W3] ) 2 

A  measure  of  quadratic  contribution  to  the  regression 
is 

nXW  *  Yxw/  *E[W41  -  1  -  e2[w3]}  Vi  . 

Note  that  ig  the  proportion  of  total  variance  of  X 

accounted  for  by  the  quadratic  component  of  the 
regression. 

To  find  the  maximum  possible  curvature,  the 
numerator  of  Yxw  is  maximized  with  respect  to  a  fixed 
b.  In  terms  of  a  and  b,  this  gives  the  non-linear 
problem  of 
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maximize  c  -  Z  arb8bt  •*  t> 

•  (Z  brbsbt  U(r,  s,  t))  •  (2  arb*  or8) 
subject  to 

§aras  a*s  *  ~  brbs  °ra  * 

Consider  c  ~  V2X  5  Ms0  rs  where  X  is  a  Lagrange 
multiplier,  and  differentiate  with  respect  to  au  to  give 
for  u  ■  1,  V  at  a  stationary  point, 

Z  b,bt  M(U,  s,  t) 

“  (  Eb^.  Z b]>bsbt  l*(f»  t) )  —  X  Zat  CTut 

■  0  . 

For  the  maximum  value  of  c  *,  they  find  that 

au  -  (I  brbs  v(r,  s,  t)  otu 

-  bu  zb^bsht  w<r,  s,  t)]/  c*  , 

where  c^j  is  the  i j^1  element  of  E  “1,  and  that  n 
the  supremum  of  over  a  for  fixed  b  is 


•  Zbfbgbtbu  u(r,s,p)  y(t,u,q)c  P<J  -  [  Zbrbsbt  U(r,s,t))3 

Zbrbsbtbu  U(r,s,t,u)  -  1  -  [E  brbsbt  y(r,a,t)]3 


This  last  expression  must  be  maximised  using  the  sample 
values  of  the  higher  moments.  The  authors  suggest  that# 
after  locating  an  initial  value  bQ,  a  "hill-climbing" 
algorithm  be  applied.  The  value  of  log  n02  is  used  as 
a  test  statistic.  This  statistic  is  tested  for 
normality. 

The  test  proposed  by  Cox  and  Small  raises 
several  questions.  Evaluating  nc2  requires  the  solu¬ 
tion  of  a  highly  non-linear  expression.  They  point  out 
the  potential  for  reaching  a  local  rather  than  global 
maximum.  They  state  that  the  effort  involved  in  evalu¬ 
ating  n o2  is  of  order  v4.  This  significantly  reduces 
the  size  of  the  matrix  Y  which  may  be  considered. 

Another  problem  for  the  user  is  that  an  inspection  of 
scatter  diagrams  for  the  derived  variables  is  required. 
This  means  that  after  nQ2  is  derived,  the  practitioner 
must  still  make  a  value  judgement.  No  indications  of 
the  power  of  this  test  are  offered. 

Small  [SMAL80]  suggested  testing  the  skewness 
and  kurtosis  of  the  marginal  distributions  of  a  multi¬ 
variate  distribution.  Suppose  the  random  variable  to  be 
tested  is  p-dimensional  and  that  a  random  sample  of 
size  n  has  been  taken  for  this  variate.  Let  Xi  be  the 
vector  of  sample  marginal  coefficients  of  skewness  wit* 


covariance  matrix  Vi,  and  Xo  be  the  vector  of  sample 
marginal  coefficients  of  kurtosis  with  covariance  matrix 
V2*  Small  then  applies  Johnson's  Su  transformation 
[  JOHN49]  component- wise  to  these  margined  sample 
statistics  to  obtain  the  vectors 

71  ■  Sisinh^fXi/  X) 

and 

Y2  *  Y2I  +  5  2Sinh"1t(X2  -  Cl)/  X2I  • 

where  S2,  Y2 »  end  C  are  found  using  the  first 

four  moments  of  Xj,  and  X2  end  tables  given  by  Johnson 
[JOHN65] .  The  components  of  yi  and  y2  have  distribu¬ 
tions  that  are  approximately  standard  normal,  so  that 
their  covariance  matrices  Ui  and  U2  have  main  diagonal 
elements  of  unity.  If  the  theoretical  correlation 
between  the  ith  and  jth  variates  of  the  original  data  is 

Pij,  the  off-diagonal  elements  of  Ui  and  02  are  given 
3  4 

by  p^  and  p^,  respectively.  The  test  statistics 

Qi  •  and  02  *  y2T?2”1af2 

each  have  an  approximate  x2p  distribution  under  the  null 
hypothesis.  Sines  the  skewness  and  kurtosis  coef¬ 
ficients  are  uncorrelated  and  nearly  independent  in  the 
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univariate  case,  Qi  and  Q2  are  nearly  independent  and 
the  significance  level  for  the  tests  may  be  determined 
accordingly.  Of  course  when  the  correlation's  {pij}  are 
unknown,  they  must  be  estimated  from  the  sample  data. 

Small's  test  offers  an  ease  of  calculation  not 
seen  in  earlier  tests.  However,  it  is  based  upon  the 
necessary,  though  not  sufficient,  condition  that  the 
marginal  distributions  must  be  univariate  normal  for  a 
multivariate  distribution  to  be  normal.  He  states  that 
is  is  possible  for  a  marked  departure  from  multivariate 
normality  to  be  accompanied  by  apparent  normality  in  the 
marginals.  Results  of  this  test  are  therefore 
inconclusive. 

Hawkins  [HAWK81]  proposed  simultaneously  testing 
for  multivariate  normality  and  homoscedasticity.  He 
employed  the  Anderson-Darling  statistic  to  test  for 
these  properties. 

Let  Xij,  i  •  1,  ...,  g,  j  ■  1,  ...,  m  be  a  vec¬ 
tor  of  p  components.  Hawkins  proposed  to  test  for 

Xij  .  Np(  5i,  Z  i> 

where  Np  indicates  a  p-variate  normal.  Let  7i,  and 
Si  respectively  denote  the  sample  mean  vector  and 
variance  matrix  of  the  sample  {Xij*  1  <  j  <  ni  >.  Let 
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N 


9 

Z  iU 
i-1 


and 


9 

S  -  l  (ni  -  l)Si/(N  -  g)  . 
i»l 


Define 


Vi;j  -  (Xij  -  Xi.)TS“l(Xij  -  Xi.)  . 

Xi«*  *nd  S*  respectively  denote  the  sample  mean 
vector  for  group  i  and  the  pooled  covariance  matrix 

obtained  by  deleting  Xij  from  the  sample,  Hawkins  shows 
that 

T2  -  (ni  -  lHXij  -  Xi.*)Ts*”1(xij  -  Ii.*)/ni 

follows  a  Ho tellings  T2  distribution  and  therefore 

Fij  -  (N  -  g  -  p)T2/[p(N  -  g  -  1)] 

has  an  P  distribution  with  p  and  N  -  g  -  p  degrees  of 
freedom.  Invoking  the  binomial  inversion  theorem 
[PRB872] ,  Hawkins  indicates  that 

Fij  -  t<N-g-p)niVij]/  {pt(ni-l)Oi-g)  -  niVijl  >  . 
The  statistic 


Aij  •  P{F  >  j ) 

is  distributed  as  a  uniform  (0/  1)  variate  under  the 
null  hypothesis.  The  following  are  the  proposed  test 
statistics.  For  all  i,  let  Ai(i)  <  Aj.(2)  <  •••  < 
Ai(nj.)  be  the  order  statistics  for  the  Aij.  The 
Anderson-Darling  test  statistic  is 


Wi  -  m  -  m-1  ^2M2j  "  l)llogAi(j)  ♦ 


log ( 1  -  Ai(ni  -  j  +  l))] 


The  Wi's  may  asymptotically  be  expressed  as 


Wi  -  I  Zik2/[k(Jc  +  1)] 
k-1 


where  Zi*  ■  -  l (2k  +  l)n]^2  Pjc  {2Ai(jj  -  1} 


and  Pfc  is  a  kth  order  orthogonal  polynomial.  TO  test 
for  nonnormality,  the  Wi  and  first  few  Sij  must  be  com¬ 
puted.  Monnormality  is  indicated  if  tho  Zij's  are  near 
zero  and  the  Wi's  are  large. 

Hawkin' s  test  is,  unfortunately,  a  qualitative, 
rather  than  quantitative  one.  There  are  no  tabulated 
critical  values  for  his  test  statistic* 
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Computationally ,  obtaining  the  Wi's  and  Ztj'a  la  a 
forbidding  task.  In  addition*  it  remains  unclear  as  to 
how  many  of  the  Z^j's  must  be  computed  and  examined. 

3.3.2  Motivation  for  the  Proposed  Test 

Shapiro  and  Wilk  [SHAP65]  concentrated  their 
efforts  on.  developing  a  test  for  univariate  normality 
which  would  be  both  scale  and  origin  invariant.  Their 
work  was  based  upon  an  attempt  to  formalize  the  depar¬ 
tures  from  statistical  linearity  of  probability  plots. 
Let  X1  £  x2  $.  *  *  *  <  *n  denote  the  order  statistics  for 
random  sample  of  size  n  from  a  standard  normal 
distribution  with 

Etxi]  »  mi 


and 


COV(Xi,  Xj)  »  Vij  . 


If  y?  »  (yi,  ...,  yn)  is  a  vector  of  order  statistics 
•  for  a  random  sample  of  sise  n  from  a  fixed  but  unknown 

*  distribution*  then  we  wish  to  test  whether  this  distri- 

bution  is  normal.  If  the  underlying  population  is 
normal*  then  we  may  write 


Let  mT  -  (mi,  mn)  and  V  *  (v^j).  For  symmetric 

distributions,  the  unbiased  least-squares  estimates  for 
p  and  a  are 

„  n 

y  *  [l/n]  Z  yi 
i-1 

and 

a  *  [mTviy] /[mTvlm]  . 

The  Shapiro-Wilk  test  statistic  is 

n  n  _ 

w  -  t  Z  *iYi)2/  Z  (yi  -  y)2 

i-1  i-1 

where 

a*  -  l«TV-l]/tmTv-lv-lm]}4  .  (3.75) 

To  compute  W  for  a  random  sample  of  size  n  *  2(1)50, 
coefficients  {a*}  are  located  in  tables  provided  by 
Shapiro  and  Wilk.  [Notes  i(j)k  denotes  the  integers  i, 
i+j,  i+2j,...,  i+»j-k.]  They  also  provide  percentage 
points  for  if  based  upon  Johnson's  [J0HN49]  St  curves. 

Shapiro  and  Wilk  compared  the  power  of  the 
W-test  with  that  of  several  standard  tests  of  normality* 
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chi-squared,  third  and  fourth  moments,  Kolmogorov- 
Smirnov,  Cramer- Von  Mises,  a  weighted  Cramer- Von  Mises 

i 

«  using  the  Anderson-Darling  method,  Durbin's  method,  and 

a  range/standard  deviation  test.  Shapiro,  Wilk,  and 

$  Chen  [SHAP68]  also  conducted  a  comparative  study  on  the 

♦ 

powers  of  several  tests  for  univariate  normality.  Both 
studies  concluded  that  the  Shapiro-Wilk  test  has  super¬ 
ior  empirical  power  over  a  wide  range  of  alternative 
nonnormal  distributions.  This  fact,  coupled  with  its 
ease  of  calculation,  mak.s  W  the  clear  choice  for  uni¬ 
variate  testing  of  normality. 

Subsequent  to  the  first  appearance  of  the 
Shapiro-Wilk  test,  two  similar  tests  have  been  sug¬ 
gested.  Tables  of  critical  values  (W  a(n)}  and  coef¬ 
ficients  (ai(n)}  for  the  Shapiro-Wilk  test  have  only 
been  provided  for  samples  up  to  size  50.  Shapiro  and 
Francia  ISHAP72]  offered  the  following  test  statistic: 

n  n 

'  W'  »  [  E  biyil2/  E  tyi  -  yl 2 

i-1  i-1 

* 

where 

* 

«  bT  •  (bj.,  b2r  ...»  bn)  ■  raT/[mTra)J£  . 


Values  for  m  have  been  determined  for  sample  sizes  of 
n  *  2(1)100(25)300(50)400,  and  the  null  distribution  of 
W'  has  been  approximated.  Weisberg  and  Bingham  [WEIS7S] 
suggested  using 

mi  =  *-l< [i  -  3/8]/ [n  +  1/4] )  , 

where  is  the  inverse  of  the  standard  normal  at 

p,  as  an  approximation  for  m  in  the  H'  statistic.  The 
W’  statistic  has  the  computational  advantage  of  no 
storage  of  constants  for  machine  use.  Both  W  and  W 
have  approximately  the  same  power  as  W  [SHAP72,  WEIS75]. 

3.3.3  A  Multivariate  Normality  Test 

Due  to  the  superior  power  of  the  Shapiro-Wilk 
test,  we  developed  a  multivariate  version  of  this  test 
based  on  the  union-intersection  principle  [MORR76]. 

Let  X  be  a  p-dimensional  column  vector  so  that 
t£i  :  1  <  i  <  n  }  forms  a  random  sample  of  size  n  from 
a  p-variate  distribution.  We  wish  to  test  the 
hypothesis  that 

H0  s  X  ~  Np(  yx,  E  x) 

for  some  mean  vector  yx  and  dispersion  matrix  £  x, 
against  the  alternative 


First  we  form  an  arbitrary  linear  combination  c^X  where 
cT  is  a  nonnull  p-diaension  real  vector.  Therefore  e®X 
is  univariate  and  we  may  test  the  univariate  hypothesis 

Hq(c)  :  c*X  -  N(ct  yx,  cT  2*c)  (3.76) 

against  its  alterative 

Hi(c)  s  cTx  -  N(cT  yx»  cT  ZxS)  . 

Let 

Yj  s  jth  order  statistic  of  {  c^X^  :  1  <,  i  ,<  n> 
J  *  *  (3.77) 

The  univariate  Shapiro-wilk  statistic  for  Y  is 

n  n  _ 

W(c)  -  [  S  SjYilVl  E  (Yj  -  X)21  •  (3.7#) 

j-1  J  3  j-1 

The  acceptance  region  for  H0(c)  is 

W(c)>Wg(n)  (3. 7#) 

where  «g  (n)  is  the  critical  value  for  the  univariate 
test.  Mow  Bo  i«  true  if  and  only  if  Self)  la  true  for 
all  nonnull  c.  Thu*  we  will  accept  s©  if  and  only  if  we 
accept  all  the  univariate  hypotheses  8©<$)  for  all  f 
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in  Rp.  Thus,  the  acceptance  region  for  Ho  is  given  by 
the  intersection 

fllW(c)  _>  W  0(n)]  .  (3.80) 

c 

This  intersection  is  equivalent  to  the  condition 

W*  *  min  W(c)  >  W  g (n)  .  (3.81) 

c 

When  the  null  hypothesis  is  true,  W*  will  have  a  par¬ 
ticular  CDF.  From  this,  critical  values  for  Wg  (p,  n) 
could  be  calculated  and  thus,  for  a  test  of  level  S  ,  we 
accept  H0  if 

W*  >  W  g (p,  n)  .  (3. 82) 

Now  locating  W*  is  equivalent  to  solving  each  of  the 
following  quadratic  programming  problems: 

QP( a  )  :  Maximize  z  *  YTY 

subject  to 

lTy  .  o 

a Ty  -  l 

*«* 

Y i  -  cTx  <J(i)  *0,  1  £  i  £  n 


where  o  (1)  is  the  image  of  i  under  the  permutation  a 


of  the  integers  1,  . ..,  n  and  a  is  given  by  equation 
(3.75).  Thus,  QP(  a  )  is  a  problem  in  n  +  p  unrestricted 
variables  with  2n  +  1  linear  constraints,  where  z  is  a 
convex  function  and  the  constraints  form  a  convex  solu¬ 
tion  space.  There  are  nl  problems  of  the  form  QP(  a  ) 
corresponding  to  each  possible  permutation  a.  We  pro¬ 
pose  the  following  heuristic  to  identify  and  solve  the 
appropriate  problem  QP(  a  ) .  Let 

n  ^ 

A  *  l  (Xj  -  XHX-j  -  X)T  .  (3.83) 

j«l  ~  ~  ~ 

Por  m  *  1,  ...,  n,  we  compute  the  statistics 

Uj  =  (Xm  -  X)TA-l(Xj  -  X)  ,  1  <  j  <  n  (3.84) 

and  let  Uq),  ...,  Un  indicate  the  ordered  values.  Mow 
we  evaluate 


W»  - - 

(Xg  —  ^)^A"^(Xg  “  7) 


1  <  m  <  n  . 


(3.85) 


Each  value  of  m  is  associated  with  a  particular 
permutation  a  defined  by  the  relation 

u(j)  ■  u  c(j)  »  1  <  j  <  n 

The  value  of  m  minimizing  (3.85)  thus  identifies  the 
permutation  aQ  and  the  problem  QP(  <j0)  which  should  be 
solved.  Ihe  final  test  statistic  is  given  by 

W*  -  (aTY0>/(*oT*o)  ■  V(*oTYo>  (3.86 

where  Y0  is  the  optimal  solution  vector  to  the 
quadrataic  program  QP(  0o).  We  must  note  that  it  has 
not  been  established  that  the  use  of  will  guarantee 
that  the  proper  permutation  is  selected.  A  rigorous 
demonstration  of  the  validity  of  this  procedure  is  the 
subject  of  ongoing  research.  Once  this  has  been 
accomplished,  the  null  distributions  of  the  resulting 
test  statistic  (3.86)  for  various  values  of  n  and  p  can 
be  estimated  by  applying  Johnson's  Sg  system  of  curves 
[JOHN49]  to  simulation-generated  empirical 
distributions. 


CHAPTER  IV 


* 


EXPERIMENTAL  PROCEDURE 

To  gain  insight  into  the  problems  of  relative 
bias  and  variance  reduction  and  to  determine  the  effec¬ 
tiveness  of  the  two-stage  control  procedure,  four 
queueing  systems  were  selected  for  simulation  testing. 
This  chapter  presents  a  description  of  these  models,  a 
discussion  of  available  performance  measures,  and  an 
analysis  of  the  validation  procedure. 

4.1  Selection  of  Experimental  Systems 

Queueing  models  have  become  a  standard  experi¬ 
mental  vehicle  for  controlled  simulation.  These  include 
open,  closed,  and  mixed  systems.  In  an  open  system, 
customers  arrive  to  the  system,  are  serviced  at  one  or 
more  stations,  and  depart.  An  example  of  an  open  system 
is  a  retail  center  where  customers  enter,  make  a 
purchase,  and  leave.  In  a  closed  system,  a  fixed  number 
of  customers  remain  within  the  system  during  its  opera¬ 
tion.  An  on-line  computer  with  a  fixed  number  of  ter¬ 
minals  is  such  a  system.  Mixed  systems  contain  both 
open  and  closed  classes  of  customers.  An  example  of 
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this  is  a  computer  system  with  batch  (open)  input  and  a 
fixed  number  of  on-line  (closed)  terminals. 

The  first  system  selected  for  analysis  was  the 
M/M/1  queue.  This  system  was  appealing  because  of  its 
analytical  tractability .  The  system  response  variables 
we  examined  were  the  customer's  total  average  time  in 
system  and  the  steady-state  proportion  of  time  that  the 
server  was  busy.  The  second  response  is  hereafter 
referred  to  as  the  station  utilization.  This  was  the 
only  open  queueing  system  we  examined. 

The  second  stochastic  model  we  chose  was  a 
periodic  review  inventory  model.  The  system  was 
operated  under  a  stationary  (s,  S)  inventory  policy. 

Let  dn  be  the  demand  in  period  n  and  Xn  be  the  total 
inventory  on  hand  at  the  beginning  of  period  n.  Thus 


Xn  +  1 


(  Xn  -  dn  ,  dn  <  Xn 
(  S  r  otherwise 


(4.1) 


where  s  is  the  reorder  point  and  S  is  the  reorder  (or 
stock  control)  level.  There  are  initially  XQ  »  s  units 
in  inventory  and  a  return  to  this  state  signals  the 
beginning  of  a  new  regenerative  cycle.  The  response 
variables  we  considered  were  X,  the  average  number  on 
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hand,  and  ir^,  the  long-run  proportion  of  time  that 
there  are  i  units  on  hand.  Examination  of  the  (s,  S) 
inventory  model  permits  us  to  study  the  behavior  of  the 
control  procedures  under  a  different  correlation 
structure.  For  a  covariance  stationary  process 
{  Xt,  t  »  1,  2,  ...}  define  pj  to  be  the  correlation 
between  X  and  X  +  j.  This  is  given  by 

Pj  -  Cov(Xt,  xt  +  j)/  °2‘ 

In  general,  the  response  time  variable  for  all  queueing 
systems  have  a  similar  correlation  structure.  Figure 
4.1  shows  the  correlation  function  for  system  1.  For 
the  (s,  S)  inventory  model  the  correlations  between  the 
inventories  on  hand  for  various  time  lags  is  negative 
for  odd-numbered  lags  and  positive  for  even-numbered 
ones.  The  correlation  function  for  this  system  is  shown 
in  Figure  4.2. 

The  third  system  we  considered  was  the  central 
server  model,  shown  in  Figure  4.3.  This  is  frequently 
used  as  a  simple  model  c  *  a  multiprogrammed  computer 
system  with  service  center  1  representing  the  processor 
and  the  other  service  centers  representing  input-output 
devices.  The  system  consists  of  three  service  centers, 
each  of  which  has  si  servers,  i  •  1,  2,  3.  The  number 
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Figure  4.4  Machine  Repair  Model 
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of  customers  (that  is*  the  level  of  multiprogramming)  is 
fixed  and  is  equal  to  N.  Initially  all  N  customers  are 
at  station  1.  Service  time  at  station  1  is  exponen¬ 
tially  distributed  with  mean  yj.  With  probability 
Pl2r  a  customer  leaving  service  center  1  immediately 
enters  service  center  2  where  he  joins  a  FIFO  queue  to 
await  service.  Service  time  at  station  2  is  iid  expo¬ 
nential  with  mean  y2.  Alternatively,  a  customer  leav¬ 
ing  station  1  enters  service  center  3  with  probability 
Pl3  *  1  -  P12  to  join  a  FIFO  queque  and  await  service 
there.  Service  times  at  station  3  are  iid  exponential 
with  mean  U3«  After  completion  of  service  at  service 
centers  2  or  3,  the  customer  returns  to  center  1.  The 
response  variable  examined  here  was  the  time  between 
successive  arrivals  by  a  customer  to  service  center  1. 

The  fourth  stochastic  model  we  examined  was  a 
variation  of  the  machine  repair  model.  This  system, 
shown  in  Figure  4.4,  consists  of  four  queues  and  a  fixed 
number  of  units  N.  Initially,  all  N  units  are  at  sta¬ 
tion  1  with  si  units  in  operation  and  N  -  si  waiting  in 
the  queue  as  spares.  The  time  until  failure  of  an 
operational  unit  at  station  1  is  iid  exponential  with 
mean  y2.  Upon  failure,  a  unit  is  sent  for  repair. 

With  probability  pi2,  a  major  repair  is  required  at 
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station  2.  With  probability  P13  ■  1  -  pi2»  a  minor 
repair  is  needed  at  station  3.  Stations  2  and  3  are 
PIPO  queues  wit i\  S2  and  33  repairmen  respectively. 

Repair  times  are  iid  exponential  with  mean  y 2  for  a 
major  repair  and  mean  y3  for  a  minor  repair.  Following 
repair,  a  unit  proceeds  to  station  4  for  inspection  on  a 
PIPO  basis  by  one  of  the  34  inspectors  each  having  iid 
exponential  service  times  with  mean  y4.  A  unit  will 
fail  inspection  with  probablity  P43.  Should  it  fail,  it 
is  sent  back  to  station  3  for  further  repair. 

Otherwise,  it  is  returned  to  station  1  where  it  joins 
the  queue  of  spares;  if  there  are  fewer  than  s^  opera¬ 
tional  units,  it  goes  into  service  immediately.  The 
response  variables  examined  were  the  average  number  of 
operational  units  at  station  1;  the  server  utilizations 
at  stations  2,  3,  and  4;  and  the  time  required  for  a 
newly-failed  unit  to  reenter  station  1. 

We  are  able  to  obtain  the  desired  values  of  the 
steady-state  parameters  analytically  for  all  four  of 
these  systems.  The  method  for  calculating  these  is 
found  in  §  4.5. 

The  first  three  models  presented  were  used  to 
analyze  the  bias  and  confidence  interval  problems  which 
exist  when  top-controlled  regenerative  analysis  is  used. 


The  fourth  model  was  selected  for  validation  of  tha  two- 
state  procedure.  Results  of  the  experimentation  appears 
in  Chapter  V. 

4 . 2  Performance  Measures 

To  determine  the  value  of  the  two-stage  method 
of  applying  concomitant  control  variables ,  some  stan¬ 
dards  of  comparison  are  required.  This  sect ion  contains 
a  discussion  of  available  performance  measures  and  a 
presentation  of  those  methods  selected  for  use. 

Most  performance  measures  compare  the  efficiency 
of  direct  simulation  (subsequently  labelled  method  0) 
to  the  efficiency  of  ease  alternative  procedure 
(subsequently  labelled  method  1)  for  variance  reduction 
or  confidence  interval  estimation.  This  relative 
efficiency  factor  of  method  1  to  method  0  may  be  based 
upon  the  coat  involved*  the  gain  or  loss  of  pceoieion 
(accuracy),  or  the  reliability  of  the  technique. 

The  cost  of  a  technique  is  generally  a 
raf lection  of  the  amount  of  computation  required  by  that 
method.  Baasersley  and  Handscomb  [HMM!  proposed  that 
this  should  be  a  measure  of  the  elapsed  coaputer  time 
t4,  i  •  o,  lr  thus  they  proposed  the 

labor  ratio  »  xi 
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as  a  measure  of  cost.  Moy  [MOY65]  suggested  using  the 
computer  processing  charges,  Y^,  required  to  obtain  a 
fixed-width  confidence  interval  by  each  technique 
i  ■  0,  1.  Hoy  used  as  his  performance  measure  the 
percentage  of  original  cost  saved, 

%  gain  -  100  (  Y0  -  Yi)/  Y0«  (4.2) 

Both  of  these  methods  are  difficult  to  use  in  practice. 
The  results  obtained  depend  on  the  computing  machinery 
used  and  the  methods  by  which  computing  time  is  charged. 

The  precision  of  method  i  is  usually  considered 
to  be  inversely  proportional  to  its  estimator  variance 
a i2.  Cochran  [COCH771  expressed  the  relative  precision 
of  method  1  to  method  0  as 

2  2 

%  relative  precision  *  100 q' al*  (4.3) 

Hammer8ley  and  Handscomb  [HAMM64]  proposed  a  comparable 
measure: 

variance  ratio  ■  o2/o2.  (4.4) 

Kish  [KXSH65]  considered  using  the  "design  effect" 


deff  -  o2/c2 


(4.5) 


This  was  later  described  by  Lavenberg  (LAVE77(I),  (II); 

78;  79]  as  the  "minimum  variance  ratio."  Kleijnen 

[KLEI74]  suggested  using  the  percentage  o£  the  variance 
2 

of  method  0,  oQ,  that  was  eliminated  using  method  1: 

%  variance  reduction  *  100 [(c|?  -  o^)/cj?)]. 

u  1  u  (4.6) 

Several  researchers  [HE1D78,  LAVE79,  1GLE79]  have 
concentrated  on  the  reduction  of  the  width  of  the 
resulting  confidence  intervals: 

confidence  interval  reduction  A  *  [tx<Ji  ]  /  [t2°2l 

(4.7) 

where  to  and  ti  are  selected  critical  values  of  the 
distributions  relevant  to  each  of  the  methods.  In  all 

A  A 

of  these  methods ,  the  sample  variances ,  oo  and  ax,  are 
used  when  necessary. 

For  an  estiraand  9 ,  the  accuracy  of  the 

A 

estimator  8  derived  by  method  i  is  an  expression  of 
the  size  of  the  error  0^—0  [COCH77,  HANS53,  RAJ68] . 
The  bias 

Bi  =  E[  8 il  -  (4.8) 

reflects  the  systematic  component  of  the  error,  while 
the  mean  square  error 
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MSEi  =  E[(§  t  -  0  )2j  «  0i2  +  Bi2  (4.9) 

takes  into  account  bias  as  well  as  precision  as  a 
measure  of  the  overall  accuracy  of  method  i.  Rao 
[RA069]  studied  the  performance  of  a  variety  of  ratio 
estimators  in  comparison  to  the  classical  ratio  estima¬ 
tor,  method  0.  He  used  two  measures  of  accuracy  for 
each  alternative  estimator: 


and 


%  relative  accuracy  for  method  i  *  lOO.HSE^/MSEo 

(4.10) 


%  bias  ratio  for  method  i 


100. 


1  Bi  1  / 


/MSEi. 

(4.11) 


Reliability  is  taken  to  be  a  measure  cf  the 
actual  coverage  probability 

Pi  *  Pr  {  0i  —  ti  ci  8  ^  0  i  +  ti  Oi  }  (4.12) 

of  confidence  intervals  for  0  that  are  constructed 
using  that  method  [WILS78;  LAVE78,  79].  Lavenberg  et 
al.  (LAVE79]  considered  the  gain  or  loss  of  coverage 
achieved  by  method  1  relative  to  method  0  under  the  same 
conditions: 

coverage  gain  *  &i  -  P&.  (4.13) 
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Wilson  [WILS78]  used  a  similar  measure  to  compose  an 
alternative  technique  but  required  that  the  confidence 
intervals  he  adjusted  to  a  common  width. 

Some  of  the  available  performance  measures 
attempt  to  combine  some  of  the  basic  measures  into  a 
single  overall  figure  of  merit.  Haramersley  and 
Handscomb  [HAMM64]  felt  that  any  measure  of  efficiency 
should  incorporate  cost  and  precision.  They  proposed 
using  a  product  of  their  variance  and  labor  ratios  to 
obtain 

efficiency  gain  -  (  tq  aJ)/(  t  x  a*).  (4.14) 

Som  [SOM73]  considered  cost  and  accuracy  to  be  of 
significance  and  offered 

relative  cost  efficiency  ■  (  Yo*mSEo)/(  Yi*MSEi) 

(4.15) 

as  a  performance  measure.  Both  of  these  standards  for 
efficiency  clearly  have  the  same  problems  associated 
with  any  cost  measure. 

The  selection  of  performance  measures  to  be 
applied  in  this  research  was  based  upon  the  decision 
that  the  evaluated  technique  should  not  disrupt  the  nor¬ 
mal  course  of  the  simulation  in  any  way.  Therefore  the 
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alternative  methods  simply  utilize  available  simulation¬ 
generated  data  in  different  manners.  Experimentation 
revealed  that  the  CPU  time  required  to  obtain  the  simu¬ 
lation  data  for  each  collection  of  experiments  greatly 

,  outweighed  the  time  involved  in  applying  the  various 

*■ 

sets  of  control  variables.  Therefore ,  any  differences 
in  cost  among  the  various  alternatives  is  overshadowed 
by  the  basic  simulation  cost.  Consequently,  cost  was 
omitted  as  a  criterion  for  the  performance  of  any 
method.  Kleijnen's  variance  reduction  percentage  (4.6) 
was  selected  because  it  seems  to  have  become  a  standard 
in  simulation  experimentation.  Although  researchers 
have  used  a  variety  of  criteria  such  as  (4.3),  (4.4), 
and  (4.5),  they  frequently  returned  to  (4.6)  in  the 
discussion  of  their  results.  Additionally,  for  control 
variate  analysis,  the  variance  reduction  percentage  cor¬ 
responds  to  the  square  of  the  coefficient  of  correlation 
[LAVE81] .  The  bias  factor  given  by  (4.8)  was  also 
•  averaged  over  several  replications  of  each  experiment  to 

determine  the  accuracy  of  the  resulting  point  estimates. 
Finally  to  incorporate  the  reliability  criterion,  the 
estimated  coverage  for  each  variant  of  controlled 
regenerative  analysis  was  computed  along  with  the 
coverage  gain  (4.13). 
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4.3  Selection  of  Experimental  Parameters 

For  each  of  the  .stochastic  systems  selected 
for  examination,  a  meta-experiment  consisting  of  50 
independent  simulation  runs  was  performed.  For  each 
experiment,  a  group  of  regenerative  tours  was  used  to 
construct  point  estimates  and  confidence  interval 
estimates  for  the  selected  response  variables.  These 
estimates  were  constructed  with  and  without  the  use  of 
control  variables.  We  then  averaged  our  selected  per¬ 
formance  measures  over  all  of  the  experiments  within  the 
meta-experiment . 

In  order  to  obtain  reliable  results  with  the 
regnerative  method,  a  fairly  large  number  of  tours  are 
required  [LAVE78] .  In  this  research,  we  used  from  500 
to  3250  regenerative  cycles  in  each  experiment.  Wilson 
[WILS79]  showed  that  when  using  top-controlled  regenera¬ 
tive  analysis,  a  large  relative  bias  is  encountered  when 
the  number  of  tours  used  is  small  (less  than  100  tours). 
Additionally,  it  has  been  found  that  a  large  number  of 
short  tours  is  preferable.  For  this  reason,  the  regen¬ 
eration  points  selected  for  each  of  our  experimental 
systems  were  those  states  which  occured  the  most  fre¬ 
quently.  In  addition,  for  the  closed  systems,  the 
number  of  customers  was  kept  lower  to  increase  the  rate 


of  regeneration.  In  order  to  obtain  the  convergence  to 
joint  normality  of  the  control  variates  for  these  short 
tours,  we  found  it  necessary  to  accumulate  controls  over 
a  batch  of  tours.  The  method  for  selecting  a  minimum 
batching  size  is  presented  in  S  5.3. 

For  the  M/M/1  queue,  light  traffic  intensity 
results  in  more  frequent  regeneration.  He  therefore 
chose  an  interarrival  rate  of  1.0  with  a  service  rate  of 
2.0  for  a  traffic  intensity  of  0.5.  The  regeneration 
epochs  are  defined  to  be  those  points  in  time  when  an 
arriving  customer  finds  the  system  empty  and  idle. 

Frequency  of  regeneration  in  the  (s,  S) 
inventory  model  is  governed  primarily  by  the  number  of 
states  and  the  demand  function.  To  limit  the  number  of 
states  we  let  s  ■  3  and  S  ■  6.  The  demand  in  each 
period  is  0,  1,  or  2  each  with  probability  1/3.  Since 
asymptotically  the  choice  of  a  regenerative  state  has  no 
effect  on  the  response  variables,  we  are  free  to  choose 
any  state  to  begin  our  cycle  [CRAN75(I)].  He  found  that 
Xn  «  3  occurs  most  frequently  and  therefore  let  X0  *  3 
and  used  it  as  our  regenerative  state. 

In  selecting  the  experimental  parameters  for 
systems  3  and  4,  the  same  criteria  were  used.  The 
regeneration  epochs  occurred  when  all  of  the 


customers/units  were  at  station  1.  Other  parameters  are 
given  in  Tables  4.1  and  4.2. 

4.4  Validation  Procedure 

Analysis  of  the  experimental  results  requires 
the  comparison  of  estimators  found  through  application 
of  the  various  methods  to  the  steady-state  values  of  the 
response  variables.  If  the  true  values  of  the  estimands 
are  known,  it  is  possible  to  estimate  the  true  coverage 
of  the  nominal  90%  confidence  intervals  derived  for 
those  parameters.  Methods  for  obtaining  the  analytical 
values  for  the  response  variables  in  each  of  the  systems 
will  be  discussed  in  this  section. 

4.4.1  Results  for  the  M/M/1  Queue 

In  a  basic  open  queueing  process,  customers 
arrive  in  accordance  with  an  interarrival  process,  join 
a  queue  to  await  service  from  one  of  the  s  servers,  and 
are  served  according  to  some  service  time  distribution. 
In  the  case  of  the  M/M/1  queue,  customers  arrive  accord¬ 
ing  to  a  Poisson  process  with  arrival  rate  X  ,  enter  a 
FIFO  queue,  and  await  service  from  one  server  whose 
service  times  are  iid  exponential  with  rate  V  .  The 
utilisation  factor  [HILL80]  for  such  a  system  is  given 
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Table  4.1 


Service  Center  Parameters  for  Systems  3  and  4 


Number 

Mean 

Number 

System 

of  units 

Service  Times 

of  Servers 

N 

»1 

U2  »3  w4 

si  s2 

S3  S4 

3 

8 

1.0 

.556  5.0 

1  1 

1  1 

4 

7 

10.0 

1.5  1.0  .5 

5  1 

1  1 

Table  4.2 

Branching  Probabilities  for  Systems  3  and  4 


PI  2 


System 


P13 


P43 


P41 


p  -  X  /  V  . 

To  obtain  the  steady-state  time  in  system , 
Little's  formula  [LITT61]  to  obtain 
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(4.15) 

employ 


W  -  L/X  -  P/C <1  -  P  )  X).  (4.16) 

Thus,  for  our  particular  system,  p  *  0.5  and  W  *  1.0. 

4.4.2  Results  for  (8,  S)  Inventory  Models 

The  results  for  (s,  S)  models  are  presented  here 
in  the  context  of  the  specific  model  selected.  System  2 
is  a  Markov  chain  with  four  states  [WAGN69] .  If  Pj.j  is 
the  probability  that  the  inventory  level  will  change 
from  i  to  j,  then  the  transition  matrix  is  given  by 


(3) 

(4) 

(5) 

(6) 

“1/3 

0 

0 

2/3" 

1/3 

1/3 

0 

1/3 

1/3 

1/3 

1/3 

0 

_  0 

1/3 

1/3 

1/3. 

(4.17) 


The  steady-state  probability  that  the  system  is  in 
state  i  is  given  by  the  solution  to  the  equations 


3  -  n 


j*  •  i 


(4.18) 

(4.19) 
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The  average  number  on  hand  ux,  la  therefore 


Z  iri.i 


(4.20) 


These  values  are  contained  in  Table  4.3. 


4.4.3  Results  for  Closed  Jackson  Queueing  Systems 

Obtaining  steady-state  values  for  the  parameters 
for  systems  3  and  4  requires  a  combination  of  computa¬ 
tional  and  theoretical  results  for  closed  Jackson  queue¬ 
ing  systems  [GORD67,  BUZE73,  SOLB78] .  Results  presented 
in  this  section  are  in  the  context  of  the  validation 
model,  system  4.  System  4  is  a  closed  Markovian  queue¬ 
ing  network  with  M  *  4  stations  and  N  ■  7  customers,  and 
lid  exponentially  distributed  service  times.  Let  P^j 
be  the  probability  that  a  unit  completing  service  at 
station  i  will  be  immediately  sent  to  station  j.  The 
routing  matrix  P  *  tPijl  is  given  by 


P  - 


0 

0 

0 

.90 


.25 

0 

0 

0 


.75 

0 

0 

.10 


0 

1 

1 

0 


(4.21) 


The  vector  of  relative  arrival  rates  X  *  [  X^,...,X  ml 
to  each  station  can  be  obtained  by  solving  the  traffic 
equation 

X  -  X  P  .  (4.22) 

Since  the  traffic  equation  only  determines  X  up  to 
scalar  multiple,  one  component  of  X  may  be  arbitrarily 
set  and  the  remaining  (H  -  1)  components  are  then 
obtained  from  (4.22).  The  relative  utilization  of 
station  i  is  given  by 

Pi  =  Xi/(si  «i),  1  <  i  <  M  (4.23) 

where  station  i  has  si  servers  each  with  service  rate 
w  i  »  1/  Ui.  The  state  space  S  of  this  system  consists 
of  the  set  of  all  M-tuples  x  «  [xi, . . . ,xm] *  where  xi  is 
the  number  of  customers  at  station  i,  1  <  i  <  H. 

Therefore 

M 

S  *  {x  :  E  xi  »  w  and  xi  is  a  non-negative 
i-1 


integer  } 


(4.24) 


*iUi> 


(4.25) 


Pi*i/[silaiXi  "  8i]  ,  xi  >  Si 


and 


G(M,N)  »  s  n  V^xi) 

xes  i«l 


(4.26) 


The  equilibrium  state  probability  distribution  for  this 
system  is  given  by  [GORD67] 

M 

ff(x)  -  [1/G(M,N)  1  n  »i(Xi)  ,  xeS  (4.27) 

i-1 

Buzen  [BUZE73]  has  developed  techniques  for  computing 
the  normalizing  constant  G(M,N) ,  the  marginal  queue 
length  distributions  for  each  station,  the  station  uti¬ 
lizations  {  Ui*,  1  <  i  <  It)  ,  and  other  performance 

measures.  The  actual  arrival  rate  X i*  to  station  1  is 
then  found  by  using  the  principle  of  job  floe  balance 

X i*  •  Oi*ai  .  (4.2ft) 

Little's  formula  (LITT411  may  then  be  applied  at  the  i*& 
station  to  obtain  tbs  steady- state  mean  waiting  time 
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Wj,*  from  the  arrival  rate  Ij.*  and  the  mean  number  of 
units  Li*  awaiting  service 

Li*  -  X i*Wi*  .  (4.29) 

Solberg's  queueing  network  analysis  program  CAN-Q 
[SOLB80]  was  used  to  compute  these  values  for  system  4. 
The  results  may  be  found  in  Table  4.4.  For  system  3  we 
are  only  interested  in  t,  the  mean  time  between 
successive  arrivals  of  a  customer  at  station  1.  The 
CAN-Q  Program  found  this  value  to  be  8.07. 

4.4.4  Results  for  the  GERT  Analysis 

For  system  4,  the  mean  response  time  is  not 
available  as  a  direct  result  of  the  CAN-Q  analysis.  To 
obtain  this  value,  we  rely  on  the  GERT  analysis  tech¬ 
niques  for  generalized  activity  networks  [PRXT66(I), 
(II),  (III);  WHIT69] . 

The  GERT  network  shown  in  Figure  4.5  represents 
the  repair  process  from  the  time  an  individual  unit  in 
system  4  fails  until  it  is  returned  to  station  1. 

Because  we  are  only  concerned  with  calculating  t,  the 
mean  response  time,  the  branch  traversal  times  are 
treated  as  constants  using  the  mean  service  times  {  14} 
found  in  Table  4.1  and  the  mean  waiting  times  {W**} 
found  in  Table  4.4.  The  branch  traversal  probabilities 


Figure  4.5  GFRT  Model  of  Response  Time  Process  for  System 
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Table  4.4 


CAN-Q  Results  for  Model  4 


Station 

i 

Utilization 

Ui* 

Mean  Queue 

Length 

Li* 

Mean  Waiting 

Time 

Wi* 

1 

4.7809 

1.004 

2.0998 

2 

.17928 

.036 

.29969 

3 

.41169 

.237 

.57678 

4 


26560 


085 


16068 


correspond  to  the  branching  probabilities  £ound  in  Table 
4.2: 


<32  “  P12 
<33  -  P13 
<34  *  P43 
<35  "  P41 


(4.30) 


Therefore  the  w-f unction  [PRIT66(I),  (II),  (III)]  for 
the  branch  from  node  j  to  node  k  is  given  by 

wjlc(y)  *  qe°y  (4.31) 


where  q  is  the  appropriate  branch  traversal  probability 
and  c  is  the  traversal  time.  Mason's  rule  [PRIT66(I)] 
holds  that  the  equivalent  w-function  for  the  open  net¬ 
work  from  node  1  to  9  is  given  by 


«E(y)  *  W12W24W46W67W78W89  +  W13W35W56W67W78W89 


1  -  W35W56W67W78W83 


(4.32) 


Combining  (4.31)  and  (4.32)  we  get 


WE(y)  ■ 

q5*«xp(  (W4*-HJ4]y)  •{<32*«xp(  lW2*+V2]y)+q3*«*P(  lW3*+V3ly>} 

1  -  q4*exp([W3*+  W3+l*4*+  M4]y) 

(4.33) 


Since  wg ( 0 )  -  1,  the  mean  response  time  is  given  by 


ug  »  3/3  y(WE(y)/WE<0)]y  ,  q 

*  q2(W2*  +  V2)  +  tq3  +  (q4/q5)HW3*  +  113]  + 

u  +  (q4/q5)HW4*  +  y4] 

«  2.5418.  (4.34) 

4.4.5  Selection  of  Ratio  Estimators. 

To  perform  our  2 -stage  procedure,  we  must  obtain 
ratio  estimators  for  each  of  the  specified  response 
variables.  A  variety  of  types  of  estimators  are 
available.  Suppose  the  response  variable  of  interest  is 
r  *  E  [f (X) ]  where  f  is  a  real-valued  function  and  X  is 
some  stationary  random  variable  associated  with  the 
simulation.  In  the  context  of  regenerative  simulation, 
we  observe  the  process  {  X(s),  s  .>  0}  ,  where  X(S)  ^  X, 
in  IID  cycles  of  lengths  {  aj  :  i  >  1}  and  we  collect 
values  {Y^  :  i  £  1  }  for  each  cycle,  where  Yi  is  given 
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(here  $i  denotes  the  ith  regeneration  epoch  so  that 
°i  *  8i  +  i  -  8i).  Under  mild  restrictions  on  the 
sample  paths  of  the  process  (X(s)  }  ,  the  regenerative 
method  ensures  that 

r  *  E[Y]/El  a )  (4.36) 

and  that  the  pairs  {  (Yj.,  0  A)  > are  iid.  Let 
UjT  m  (Yi#Ci)  be  a  column  vector  with  mean  vector  y  and 
covariance  matrix  I.  Given  a  set  of  n  cycles,  we 
denote  the  sample  mean  vector  by  U  and  the  sample 
covariance  matrix  by  £  ■  (  3ij].  The  classical 
regenerative  ratio  estimator  is  given  by 

rc(n)  -  Y/a  (4.37) 

Fieller  [FIEL40]  proposed  using  the  estimator 
Y  a  -  Jc  c12 

rf(n)  -  -  *  (4.38) 

a  2  -  k  $22 

where  k  »  z2q  _  Y/2)/n«  This  estimator  is  the  midpoint 
of  Fieller' s  100  (1  -Y  )%  confidence  interval.  Three 
other  ratio  estimators  have  been  constructed  in  an 
attempt  to  reduce  the  bias  found  in  the  classical 
method.  The  jacknife  estimator,  an  extension  of  the 
work  of  Quenouille  (QUEN49,  56]  is  given  by 
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* 


4 


11 

r,(n)  -  [1/n]  Z  [n(Y/o  )  -  (n-l)(  Z  Y*/  Z  at)]. 
J  i-1  k*i  kj»i 

(4.39) 

Tin  (TZN60]  proposed 

rt(n)  *  (Y/o  ]  (l+[  Si2/(Y  a)  -  $22/  a2]/n}  . 

(4.40) 

Beale  [BEAL62]  offered  a  similar  estimator 

rb(n)*(Y/  a].{[l+  312/(nY  a) ] / [1  +  o22/(n  a2)]  }. 

(4.41) 

All  of  these  ratios  are  strongly  consistent  and  are 
biased.  Iglehart  [IGLE75]  compared  these  estimators  in 
their  use  in  regenerative  simulation.  For  long  simula¬ 
tion  runs  (i.e.  many  tours,  as  we  have  used  here),  he 
concluded  that  the  classical  ratio  estimator  is  the  pre¬ 
ferred  choice.  We  therefore  chose  to  use  that  type  of 
estimator  in  this  research. 

We  shall  now  proceed  to  show  the  specific  ratios 
employed.  For  the  (s,  S)  inventory  model  we  define  the 
following  functions: 

fl(i)  -  i 

(4.42) 

f  A  (  i )  m  ^(4)(^)  t  ^  ■  3,... ,6 
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where  l(i)  is  the  indicator  function  for  state  A.  It 
$i  is  the  starting  tine  for  the  ith  regenerative  cycle, 
and  ai«&i+i-®i,  then  we  have  for  the  “reward" 
on  the  ith  cycle 

ei+l  -  1 

Yi<*)  -  Eff£(Xj)  ,  i  >  1,  A-  1,  3,  4,  5,  6. 

j  -  Bi  “  (4.43) 

Crane  and  Iglehart  [CRAN74]  showed  that 

E  [  f  X )  ]  -  E(YiU)]/B[  cijJ  .  (4.44) 

Therefore,  our  ratio  estimators  are 

X  «  Y(l)/a 

,  (4.45) 

»  Y(A)/a  ,  A»  3,.. .,6 

For  system  3,  the  steady-state  parameter  of 
interest  is  RT*,  the  mean  response  time.  Let  t£  denote 
the  ith  observation  of  response  time — that  is,  the 
ith  time  between  successive  arrivals  by  a  customer  to 
station  1.  Let  Y£  •  n  •  (  »i)  where  N  is  the  fixed 
number  of  customers  in  system  and  a£  is  the  ith  cycle 
length.  Define 

Ai  •  number  of  arrivals  to  station  1 
on  the  i^h  tour. 


From  Little's  formula  and  the  regenerative  structure  of 
the  system  it  can  be  shown  that  [LAVE77c] 


n 

RT*  «  lim[l/n]  Z  ti  w.p.  1 
n-**»  i«l 

-  E[Yil/E[Ai]  . 

This  sugests  that  we  use 

RT  *  IN  •  a] /A 


(4.46) 


(4.47) 


as  our  regenerataive  estimator,  where  a  and  A  are  the 
sample  means. 

For  system  4,  the  steady-state  parameters  to  be 
estimated  are  the  average  response  time  RT*  from 
failure  until  completion  of  repair  and  the  station  uti¬ 
lizations  {  Oi*,  1  <,  i  <,  4  }  .  Let  xi(t)  be  the  number 
of  customers  at  station  i  at  time  t.  We  define  the 
following  variables  on  the  kth  regenerataive  cycle: 

Rjc  -  number  of  units  completed  on  the  kth  tour 

Gjc(i)  ■  /  min  {  x*(  t) ,  si  >  dt,  1  <  i  <  M 
0k  ~ 

*  0k+l 

H*  ■  Z  /  Xi(t)dt  . 

i-2  0* 
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Prom  Little's  formula  and  the  regenerative  structure  of 
the  system  it  can  be  shown  that  [LAVE77c] 


n 

RT*  *  limll/n]  I  ti  w.p.  1 

n-H»  i»l 


-  ElYil/EtAiJ  . 


(4.46) 


This  sugests  that  we  use 

RT  -  [N  •  a] /A  (4.47) 

as  our  regenerataive  estimator*  where  a  and  A  are  the 
sample  means. 

For  system  4,  the  steady-state  parameters  to  be 
estimated  are  the  average  response  time  RT*  from 
failure  until  completion  of  repair  and  the  station  uti¬ 
lizations  {  Uj.*,  1  <  i  <  4  }  .  Let  xi(t)  be  the  number 
of  customers  at  station  i  at  time  t.  We  define  the 
following  variables  on  the  kth  regenerataive  cycle: 

R]c  ■  number  of  units  completed  on  the  kth  tour 

Gjc(i)  »  /  min  (  xi(t) ,  si  >dt,  1  <  i  <  M 
8k  ”  “ 

M  6*+! 

H*  »  Z  /  Xl(t)dt  . 

i-2  9* 
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From  Little's  formula  and  the  regenerative  structure  of 
the  system  it  can  be  shown  that  [LAVE77c] 


RT*  •  lim[l/n]  I  ti  w.p.  1 
n-N»  i»l 


E[Yil/E[Ail 


(4.46) 


This  sugests  that  we  use 


RT  ■  [N  •  aJ/A 


(4.47) 


as  our  regenerataive  estimator,  where  a  and  A  are  the 
sample  means. 

For  system  4,  the  steady-state  parameters  to  be 
estimated  are  the  average  response  time  RT*  from 
failure  until  completion  of  repair  and  the  station  uti¬ 
lizations  {  Ui*,  1  £  i  <  4  }  .  Let  xi(t)  be  the  number 
of  customers  at  station  i  at  time  t.  We  define  the 
following  variables  on  the  kth  regenerataive  cycle: 

R]t  ■  number  of  units  completed  on  the  kth  tour 
,  8jc+1 

Gfc(i)  -  /  min  {  xi(t) ,  si  >  dt,  1  <  i  <  M 
8k  " 

*  8k+l 

Hfc  »  Z  /  xi(t)dt  . 
i-2  Bfc 
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Thus  Gk(i)  is  the  time-integrated  number  of  busy  servers 
at  station  i  during  cycle  k  and  Hk  is  the  time- 
integrated  number  of  units  being  repaired  in  the  kth 
cycle.  Using  the  standard  regenerative  argument,  we 
find 

ui*  *  E[Gi(i)]/E[  oil  »  1  <  i  <  M  .  (4.48) 

which  is  estimated  by 

A  __ 

ui  -  G(i)/a,  1  <  i  <  M  .  (4.49) 

To  estimate  RT*  we  again  employ  Little's  formula. 
Consider 


E[Hkl/E(akl  ■  expected  number  of  customers 
undergoing  repairs 


steady-state 
completion  rat 


expected 
response  time 


E[Rk) 

-  .  RT*  . 

Elojc) 


Thus  we  have 

RT*  -  E [Hk]/E(Rk]  ,  (4.50) 

which  implies  that  our  regenerative  estimator  should  be 
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RT  -  H/R  ,  (4.51) 

where  H  and  R  indicate  the  sample  means  over  a  set  of 
regenerative  tours. 

Simulation  models  of  the  various  queueing 
systems  were  coded  in  SIAM  [PRIT79].  FORTRAN  IV 
routines  were  used  to  perform  the  control  variable 
procedures.  Listings  for  these  programs  appear  in  the 
Appendix.  The  results  of  the  experimentation  and  an 
analysis  of  their  implications  are  found  in  Chapter  V. 


CHAPTER  V 


EXPERIMENTAL  RESULTS 

This  chapter  presents  the  results  of  the 
meta-experiments  described  in  Chapter  IV.  The  first 
section  contains  the  results  of  the  first  three  models 
using  top-controlled  or  bottom-controlled  regenerative 
analysis.  The  results  are  discussed  and  their  implica¬ 
tions  for  the  two-stage  method  are  presented.  In  the 
second  section  we  present  the  results  of  the  fourth 
system,  our  validation  model,  and  we  compare  the  two- 
stage  procedure  with  the  other  techniques.  The  last 
section  summarizes  the  findings  of  this  research  and 
presents  guidelines  for  the  practical  application  of  the 
method. 

5.1  Demonstration  of  Bias  and  Coverage  Problems 

Throughout  the  literature  we  have  seen  a 
multitude  of  attempts  to  control  the  numerators  of 
various  ratio  estimators.  Iglehart  and  Lewis  [IGLE79] 
mentioned  the  possibility  of  controlling  the  denominator 
but  did  not  pursue  the  idea.  In  this  section  we  present 
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the  results  of  systems  1,  2,  and  3  where  we  attempted  to 
explore  the  merits  of  both  of  these  techniques. 

To  effectively  control  the  numerator  or 
denominator  of  a  ratio  estimator,  we  are  faced  with  task 
of  finding  controls  which  are  strongly  correlated  with 
the  appropriate  regenerative  measurements.  The  use  of 
the  "standardized  service-time"  and  "standardized  flow" 
variates  appears  to  fill  this  requirement  while  taking 
advantage  of  all  of  sampling  performed  during  the  normal 
course  of  a  simulation. 

In  Wilson's  [WILS79]  work  he  found  that  it  was 
necessary  to  group  the  regenerative  tours  into  batches 
to  insure  the  convergence  to  normality  for  the  service¬ 
time  variables.  In  this  research  we  found  that  batching 
was  also  necessary  to  obtain  the  convergence  to  joint 
normality  of  both  classes  of  controls.  Applying  the 
univariate  Shapiro-Wilk  test  to  each  control  separately 
in  an  overall  Bonferroni-type  test  for  joint  normality, 
we  determined  a  minimum  batching  size  which  was  used 
throughout  the  experimentation.  This  batching  procedure 
is  fully  explained  in  $  5.3.1.  Table  5.1  displays  the 
batch  sizes  used.  The  observations  {(Ykrxk>  *  1  £  k  £  n> 
for  each  tour  were  averaged  over  the  n/v  batches  of 
size  v: 
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Table  5.1 


Tour-Batching  Used  in  Experimentation 


System 


Tours/Batch 


Batch/Experiraent 


jv 

Yi(v)  -  [1/v]  •  £  Yi 

J  i»(j-l)v+l 

1  <  i  <  n/v  *  k  (5.1) 
jv 

Xj(v)  -  [1/v]  •  £  Xi 

J  i-(j-l)v+l 

For  the  jth  batch  of  v  tours,  the  standardized  controls 

A.  a  [Cj,  Dj]  (5.2) 

were  collected  over  the  QC  work  stations  and  QD 
branching  points.  In  the  case  of  top-controlled 
analysis,  we  performed  a  regression  of 

Zj(v)  -  Yj { v)  -  rXj(v)  (5.3) 

on  the  components  of  &  for  each  of  the  response 
variables.  For  bottom-controlled  analysis,  we  performed 
a  regression  of  Xj(v)  on  the  components  of  $.  For  each 
of  the  first  three  systems  we  shall  use  the  following 
labelling  scheme:  (1)  meta-experiments  labelled  "A” 
refer  to  uncontrolled  estimation;  (2)  "B"  meta¬ 
experiments  are  those  using  top-controlled  estimators; 
and  (3)  "C"  meta-experiments  are  those  using  bottom- 
controlled  estimators. 

For  system  1,  we  have  two  sampling  procedures 
involved:  the  sampling  of  interarrival  times  and  ser¬ 
vice  times.  For  each  of  these,  service-time  controls 
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were  constructed  and  applied  to  the  numerator  and  the 
denominator.  Tables  5.2,  5.3,  and  5.4  present  the 
,  results  o£  our  selected  performance  measures.  We  see  a 

marked  bias  problem  in  our  top-controlled  estimates.  We 
»  also  find  large  variance  reductions  accompanied  by  great 

4 

degradations  in  coverage.  In  the  case  of  bottom- 
controlled  analysis,  we  find  that  we  have  over-corrected 
the  bias  in  the  classical  estimate.  While  the  con¬ 
fidence  interval  coverage  is  acceptable,  it  may  be  a 
direct  result  of  the  variance  increase. 

In  the  (s,S)  inventory  model,  the  only  sampling 
which  occurs  during  this  simulation  is  that  of  the 
demand  distribution.  Thus  our  control  variable  for 
system  2  is  based  upon  the  periodic  demand  d£.  Tables 
5.5,  5.6,  and  5.7  display  the  values  of  the  performance 
measures  for  this  system.  Again,  we  find  that  the  top- 
controlled  estimator  is  generally  biased.  The  control 
has  very  little  correlation  with  (Zj  :  1  <.  j  £  n/v>. 

#  Hence,  we  find  very  little  variance  reduction  in  meta¬ 

experiment  B  and  consequently  there  are  no  problems  in 
confidence  interval  coverage.  Applying  the  control 

» 

variable  to  the  denominator  produced  somewhat  mixed 
results.  In  general,  the  bias  problem  is  worsened  by 
the  use  of  the  bottom  control,  although  in  the  case  of 
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Table  5.2 


Bias  in 

Ratio  Estimators 

for  System  1 

Meta-Experiment 

Estimand 

RT 

U 

A 

.0052 

.0015 

B 

.0536 

.0148 

C 

-  .0058 

-  .0014 

Table  5.3 


Variance  Reduction 

Percentages 

Achieved  in  System  1 

* 

Me ta-Exper iment 

RT 

Estimand 

U 

♦ 

B 

57 

95 

* 

C 

-  33* 

-  132* 

♦ 

♦indicates  a  variance 

increase 

! 


i 
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Table  5.4 

Coverage  of  Nominal  90%  Confidence 
Intervals  for  System  1 


% 

4 

Meta-Experiment 

RT 

Estimand 

U 

A 

89 

88 

B 

28* 

0* 

f 

f 

C 

84 

88 

♦Significantly  below  the  90%  level 


Table  5.5 

Bias  in  the  Point  Estimator  for  System  2 


Me ta-Exper imen t 

X 

*3 

Estimand 

*4 

*5 

*6 

A 

-  .0012  - 

.0002  - 

.0004 

.0026  - 

.0020 

B 

.0019 

.0001  - 

.0009 

.0033  - 

.0025 

C 

-  .0230  - 

.0012  - 

.0017 

.0018  - 

.0036 

Table  5.6 


Variance  Reduction  Percentages  Achieved 
in  System  2 


Meta-Experiment 

X 

*3 

Estimand 

*4 

*5 

B 

3 

2 

3 

7 

7 

C 

-  139* 

-  23* 

-  10* 

25* 

-  19* 

^Indicates  variance  increase 


Table  5.7 


Coverage  of  Nominal  90%  Confidence 
Intervals  for  System  2 


Meta-Experiment 


Estimand 
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*5  there  was  a  marked  improvement.  As  was  observed  in 
system  1,  valid  confidence  intervals  were  obtained  under 
the  condition  of  variance  increases. 

After  examining  the  results  of  the  first  two 
systems ,  we  determined  that  due  to  their  lack  of  corre¬ 
lation  with  the  denominator  alone,  the  standardized 
service-time  variables  have  very  little  effect  when  used 
in  the  denominator.  Moreover,  in  a  broad  range  of 
queueing  systems,  these  variables  have  demonstrated 
their  value  as  top  controls.  We  therefore  chose  to  con¬ 
sider  service-time  variates  for  controlling  the 
numerator  and  flow  variates  for  use  in  the  denominator. 

In  system  3  we  have  available  three  top  controls 
(service-time  variables)  and  one  bottom  control  (flow 
variable) .  Results  for  the  selected  performance  mea¬ 
sures  appear  in  Tables  5.8,  5.9,  and  5.10,  and  are 
similar  to  those  found  in  the  other  systems. 

Top-controlled  regenerative  analysis  appears  to 
be  plagued  with  two  fundamental  defects.  First,  the 
procedure  causes  a  relatively  large  bias  to  be  intro¬ 
duced  into  the  ratio  estimator.  Second,  the  variance 
reductions  obtained  are  over-estimated  sinee  the  vari¬ 
ance  estimator  appears  to  systematically  underestimate 
the  true  variance.  Taken  together,  these  phenoaama 


Table  5.8 


Bias  in  the  Point  Estimator  for  System  3 


Meta-Experiment 

Estimand 

RT 

A 

.0025 

B 

-  .0685 

C 

-  .2931 

Table  5.9 


Variance  Reduction  Percentages  Achieved 
in  System  3 


Meta-Experiment 


Estiaand 

RT 


B  82 

C  -  371* 


* Indicates  variance  increase 
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Table  5.10 

Coverage  of  Nominal  90%  Confidence 
Intervale  for  System  3 


Meta-Experiment 

Estimand 

« 

RT 

!• 

A 

94 

B 

40* 

C 

58* 

•Significantly  below  the  90%  level 


» 
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result  in  degradation  of  confidence  interval  coverage. 
Schruben  [SCHR79]  evaluated  several  causes  of  loss  of 
coverage  in  confidence  intervals.  If  a  normally 
distributed  estimator  has  bias  B  and  variance  2,  he 
showed  that  the  coverage  of  a  nominal  100(1  -  )% 

confidence  interval  is  given  by 

*(-B/o  +  2X_  d/2)  -  ♦  ( -B/ a  -  2i-a/2)  <5.4) 

where  ♦  (•)  is  the  standard  normal  distribution 
function.  He  also  found  that  a  bias  as  small  as  one 
tenth  of  the  standard  deviation  creates  coverage 
problems.  In  addition,  Schruben  showed  that  the  perfor¬ 
mance  of  the  interval  estimator  decreases  as  the  sample 
size  is  increased.  In  terms  of  regenerative  simulation, 
he  found  that  when  the  point  estimate  contains  little 
bias,  the  confidence  intervals  tend  to  be  too  wide;  when 
a  large  bias  exists,  the  interval  widths  are  too  narrow. 
Thus,  we  see  that  if  the  bias  is  increased  and  we  obtain 
a  simultaneous  variance  reduction,  the  coverage  will 
fall  significantly  below  the  nominal  level. 

In  the  context  of  bottom-controlled  estimators, 
other  problems  are  in  evidence.  The  application  of  bot¬ 
tom  controls  appears  to  over-correct  the  bias  in  the 
classical  estimator.  This  may  in  turn  lead  to 
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increasing  the  magnitude  of  the  bias.  When  this  occurs, 
the  coverage  declines  as  in  the  case  of  top-controlled 
estimators.  A  greater  problem  in  bottom-controlled 
estimators  is  the  accompanying  increases  in  the 
variance.  While  this  is  generally  undesirable,  we  see 
from  equation  (5.4)  that  such  an  increase  will  maintain 
the  prescribed  coverage  provided  the  bias  is  not 
significantly  increased. 

5.2  Experimental  Results  for  the  Two-Stage  Procedure 

The  variance  reductions  achieved  using  the  top- 
controlled  regenerative  estimators  for  systems  1  and  3 
ranged  from  57%  to  95%.  If  the  top-controlled  bias 
problem  could  be  solved  and  the  confidence  interval  cov¬ 
erage  improved,  the  technique  would  yield  practical  and 
beneficial  results.  Wilson  [WILS79]  suggested  a  poten¬ 
tial  solution  to  these  problems.  By  increasing  the 
batch  size  while  holding  the  total  number  of  batches 
constant,  he  found  that  coverage  could  be  significantly 
improved.  This  is  confirmed  by  the  development  in 
Chapter  III  which  shows  that  the  bias  in  top-controlled 
point  estimators  is  of  order  1/  /n.  This  technique, 
however,  contradicts  the  purpose  of  applying  controls) 
we  wish  to  gain  more  information  from  shorter  simulation 
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runs  rather  than  being  forced  into  longer  and  more 
expensive  simulations. 

When  the  experimental  results  of  the  top-  and 
bottom-controlled  estimators  are  examined  together,  we 
see  that  the  strengths  of  one  technique  are  the  weak¬ 
nesses  of  the  other.  Thus,  taken  together,  the  methods 
may  be  able  to  compensate  for  each  other's  deficiencies. 
The  implementation  of  this  idea  is  the  two-stage 
estimator. 

In  the  validation  model,  system  4,  there  are 
four  queues  which  give  rise  to  four  service-time  varia¬ 
bles,  and  two  branching  points  which  enable  us  to  form 
two  standardized  flow  variables.  For  the  two-stage 
procedure,  we  previously  decided  to  apply  all  of  the 
service-time  variables  to  the  numerator;  for  the  denom¬ 
inator,  the  flow  variate  with  the  larger  correlation  was 
selected.  The  observed  results  {(Yi,  Xj.)  :  1  <  i  <  nJ 
of  the  n  ■  1200  simulated  tours  for  each  experiment  were 
averaged  over  k  *  50  batches  each  of  size  v  -  24.  For 
the  jth  batch  of  tours,  we  accumulated  the  standardized 
components  of  the  corresponding  control  vectors 

Cj  -  ICij,  C2j,  ...»  CQc,jlT  (5.5) 

DQD,jJT 


Pj  -  (Dlj,  D2j,  ...» 


(5.6) 
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where  QC  and  QD  are  the  number  of  top  and  bottom 
controls,  respectively.  Next  we  selected  the  column  of 
the  matrix  D  which  had  the  strongest  correlation  with 
X(v),  say  D*,  and  performed  a  regression  of  Xj(v)  on 
Dj*.  Let 

Xj*<v)  «  Xj ( v)  -  dDj*  (5.7) 

where  d  is  the  regression  coefficient.  We  next 
performed  a  regression  analysis  of 

Zj  -  Yj(v)  -  rX j ( v )  (5.8) 

on  the  components  of  Cj,  1  <  j  <  k.  This  procedure 
was  performed  for  each  of  the  response  variables  of 
interest.  In  addition,  we  also  determined  the  uncon¬ 
trolled,  top-controlled  and  bottom-controlled  estimators 
and  the  corresponding  confidence  intervals  to  allow  a 
complete  comparison  of  the  three  techniques  for  con¬ 
trolled  regenerative  analysis.  Table  5.11  indicates 
which  controls  were  applied  in  each  run.  Controls  1,  2, 
3,  and  4  are  the  service-time  variates,  and  controls  5 
and  6  are  the  flow  variates.  The  correlations  of  con¬ 
trols  5  and  6  with  the  number  completed  per  cycle  were 
-  .046  and  .063,  respectively.  Their  correlations  with 
the  tour  length  were  -  .114  and  .047  respectively. 
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Thus,  for  the  two-stage  estimator,  control  6  was  applied 
to  the  denominator  of  the  RT  estimator  and  control  5  was 
used  for  the  utilization  estimators. 

Tables  5.12,  5.13,  and  5.14  present  the  results 
of  the  meta-experiments.  The  biases  observed  in  the 
uncontrolled  and  one-stage  controlled  estimators  closely 
resemt  e  those  found  in  systems  1,  2,  and  3.  In  every 
case  but  Uj.,  the  bias  found  in  the  two-stage  estimator 
is  smaller  than  that  of  the  top-controlled  estimate. 

(An  explanation  of  the  problems  with  Ui  will  be  pre¬ 
sented  in  the  next  section. )  The  observed  variance 
reductions  and  increases  for  meta-experiments  B  through 
E  are  also  similar  to  those  found  earlier.  Applying  the 
two-stage  method  appears  to  yield  variance  reductions 
which  are  similar  although  slightly  smaller  than  those 
found  in  top-controlled  analysis.  Prom  a  decision¬ 
maker's  viewpoint,  perhaps  the  most  important  aspect  of 
a  simulation  is  the  confidence  interval.  It  is  here 
that  the  two-stage  method  proves  its  merit.  For  every 
estimand,  the  technique  raises  the  coverage  of  the  top- 
controlled  estimator,  giving  truly  valid  confidence 
intervals. 

To  summarize  the  findings  thus  far  in  this 
chapter,  we  have  seen  that  top-controlled  regenerative 
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Table  5.11 


Standardized  Control  Variates  Selected 
for  Use  in  System  4 


Meta-Experiment 

Selected  Controls  Type  of  Estimator 

A 

none 

classical 

B 

1#  2 

»  3,  4 

top-controlled 

C 

5 

bottom-controlled 

D 

6 

bottom-controlled 

E 

5 

,  6 

bottom-controlled 

F  1, 

2,  3,  4 

and  (5 

or  6)  two-stage 

Table  5.12 

Bias  in  the  Ratio 

Estimators  for  System  4 

Meta-Experiment 

Estimand 

RT 

Ul 

U2  U3  U4 

A 

.0025 

-  .0004 

-  .0010  .0006  -  .0004 

B 

.0158 

-  .0045 

.0005  .0034  .0014 

C 

-  .0011 

-  .0022 

-  .0011  .0005  -  .0006 

D 

-  .0006 

-  .0062 

-  .0013  .0001  -  .0008 

E 

-  .0048 

-  .0090 

-  .0014  -  .0001  -  .0009 

F 

.0126 

-  .0103 

.0003  .0029  .0011 
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Table  5.13 


Variance  Reduction  Percentages  Achieved  in  System  4 


Meta-Experiment 

RT 

°1 

Estimand 

U2 

1*3 

U4 

B 

64 

52 

61 

78 

89 

C 

-  9* 

-  136* 

-  8* 

2 

-  5* 

D 

-  23* 

-  283* 

-  2* 

16* 

-  18* 

E 

-  31* 

-  409* 

-  9* 

14* 

-  23* 

F 

59 

-  214* 

60 

62 

73 

•indicates  variance  increase 


Table  5.14 

Coverage  o£  Nominal  90%  Confidence 
Intervals  for  System  4 


Meta-Experiment  Estimand 


RT 

Ul 

U2 

«3 

U4 

A 

92 

92 

92 

97 

92 

B 

84 

82 

82 

80* 

82 

C 

90 

88 

96 

86 

92 

D 

90 

96 

94 

88 

92 

E 

94 

96 

98 

86 

94 

P 

90 

90 

92 

88 

88 

•Significantly  below  the  90%  level 


estimators  provide  significant  variance  reductions 
but  ultimately  result  in  unsatisfactory  coverage  of 
confidence  intervals.  The  bottom-controlled  method  will 
provide  valid  intervals,  but  results  in  variance 
increases.  The  two-stage  technique  has  been  able  to 
achieve  large  variance  reductions  while  maintaining  the 
nominal  level  of  coverage. 

5.3  Guidelines  for  Using  the  Developed  Procedures 
This  section  consolidates  the  findings  of 
sections  5.1  and  5.2  into  practical  guidelines  for  the 
use  of  the  developed  technique.  Test  procedures  are 
also  presented  to  aid  the  practitioner  in  avoiding 
potential  problems. 

5.3.1  Insuring  Convergence  of  Concomitant  Var^atvs 
The  two-stage  method  is  dependent  upon  the 
convergence  of  the  control  variates  to  joint  normality. 
In  practice,  we  must  insure  that,  the  vector  of  controls 
is  sufficiently  close  to  the  limiting  multivariate  nor¬ 
mal  distribution  for  the  expressions  of  relative  bias 
and  variance  to  be  valid.  To  insure  adequate  con¬ 
vergence  of  the  concomitant  variables,  the  variates  must 
be  accumulated  over  time  periods  long  enough  for  the 
sample  sizes  observed  at  the  work  stations  and  branching 
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points  to  produce  a  central-limit  effect.  While  we  have 
proposed  a  multivariate  Shapiro-Willc  test  in  Chapter 
III,  tables  of  critical  values  are  not  yet  available. 
Thus,  we  are  as  yet  unable  to  apply  this  procedure.  As 
an  alternative,  we  recommend  the  following  procedure 
based  upon  the  Bonferroni  inequality  and  the  Wiesberg 
and  Bingham  [WEIS75]  version  of  the  Shapiro-Wilk  test: 

1.  Based  upon  cost  and  feasibility,  choose  a  sample 
size  n  representing  the  number  of  cycles  to  be  simulated 
in  a  pilot  run.  If  possible,  take  n  1000. 

2.  For  each  of  the  n  cycles  accumulate:  (1)  nij, 
the  number  of  service  times  started  at  the  ith  service 
center  or  customers  passing  through  the  (i  -  QC)th 
branch  during  cycle  j;  and  (2)  the  raw  (unstandardized) 
controls 


CRi-;  «  E  Pijk  ,  1  <  i  <  QC  +  QD,  1  <  j  <  n  (5.9) 
J  k-1  J  -  -  -  - 


where  Pijk  is  the  kth  sample  drawn  at  control  point  i  in 
cycle  j. 

3.  Let  v  be  the  batching  size,  t  the  test  variable 
index,  and  i  the  number  of  batches.  Set  v  ■  t  ■  1  and 
l  ■  n. 
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4.  Compute  the  ^dimensional  vector  m  whose  jth 
component  is  given  by 


mj  -  4-l<  (j  -  3/8}  /  U  +  1/4}  )  ,  1  <  j  <  ft  . 

“  T5.10) 


5.  Compute  the  standardized  controls 


Cijy  -  (nijv)'^*  Z  (CRik  -  nijcu-)/  a.  (5.H) 
k»(j-l)v+l 


where 


jv 

njiv  -  _  S  «ilc  »  1  <  j  <  4  ,  (5.12) 

k«(j-l)n+l 

and  let  C^v  denote  the  sample  mean  for  the  ith 
standardized  control  using  batch  size  v. 

6.  Compute  the  modified  Shapiro-Wilk  statistic  for 
the  tth  control 

( ( t ) v ) 2/ ( ) 

w'  *  “T - —  »  (5.13) 

j«l<Ctjv  “  7tv>2 

where  C(t)V  indicates  the  sorted  statistics. 

7.  For  a  given  significance  level  a,  compare  W*  to 
the  critical  value  W'(  o/Q) ,  Q  »  QC  +  QD,  given  in 
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(SHAP72] .  If  W  >W'(  a/Q),  then  variable  t  may  be 
regarded  as  univariate  normal.  Set  t  ■  t  +  1.  If 
t  >  Q,  the  procedure  terminates/  otherwise/  repeat 
step  6.  If  W  <WM  o/Q),  let  v  -  v  +  1  and  1-  [n/v] , 
where  [a]  is  the  greatest  integer  function/  and  return 
to  step  4. 

This  procedure  will  insure  that  we  have  selected 
a  batch  size  v  large  enough  for  the  control  variables  to 
be  simultaneously  univariate  normal.  With  a  batch  size 
determined,  the  total  number  of  cycles  required  for  each 
experiment  may  be  evaluated  using  the  formula  N  -  v*nv, 
where  nv  is  the  number  of  batched  observations  desired. 

5.3.2  Selection  of  Controls 

The  objective  in  applying  the  two-stage 
procedure  is  to  obtain  a  variance  reduction  for  the 
estimator  while  keeping  bias  to  a  minimum.  Equations 
(3.72)  and  (3.74)  for  the  relative  bias  and  variance  of 
the  ratio  estimator  r(b,  d)  form  the  basis  for  our 
recommendations. 

Due  to  the  loss  factors,  it  is  recommended  that 
the  numbers  of  controls  QC  and  QD  be  kept  low.  For  QD, 
we  suggest  that  one  flow  variable  be  selected  for  the 
denominator.  The  flow  variate  which  h*s  the  largest 
correlation  with  the  denominator  is  the  logical  choice. 
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To  select  variables  to  apply  to  the  numerator  of  the 
regenerative  ratio,  Wilson  [WXLS79]  gives  a  procedure 
based  upon  step-wise  regression.  We  will  point  out  that 
this  is  also  based  upon  choosing  a  vector  which  will 
result  in  the  largest  correlation  with  the  numerator. 

5.3.3  Follow-Up  Analysis 

After  the  regression  procedures  have  been 
performed,  some  subsequent  analysis  is  required.  In 
system  4,  we  saw  that  we  had  achieved  our  stated 
objective  with  each  of  our  estimands  but  Ui.  For  the 
practitioner  it  is  of  vital  importance  to  evaluate  the 
magnitude  of  the  relative  bias  and  to  examine  the 
variance  of  the  two-stage  estimator.  If  he  has  followed 
the  recommendations  of  S  5.3.1  and  5.3.2,  he  has  only  to 
examine  two  correlations  to  determine  if  the  two-stage 
procedure  has  resulted  in  improved  performance  of  the 
ratio  estimator. 

To  evaluate  the  new  variance  estimator,  the 
correlation  between  Y*,  the  controlled  numerator,  and 
X*,  the  controlled  denominator,  must  be  examined.  If 
this  correlation  is  positive,  a  variance  reduction  will 
be  realized.  Should  this  factor  be  negative,  however,  a 
variance  increase  is  likely.  The  magnitude  of  a  nega¬ 
tive  correlation  will  influence  the  size  of  the  increase 
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(i.e.,  small  correlations  give  little  variance 
increase).  In  the  case  of  Ui  in  system  4,  Y*  and  X* 

Show  a  strong  negative  correlation.  This  may  be  seen 
from  the  fact  that  fewer  busy  servers  at  station  1  imply 
that  more  units  are  under  repair*  thus  causing  an 
increase  in  cycle  length.  This  correlation  between  Y* 
and  X*  gives  rise  to  the  variance  increase  for  Ui. 

After  examining  the  correlation  between  Y*  and 
X**  the  practitioner  should  next  evaluate  the  relative 
bias.  Again,  if  the  previous  suggestions  have  been 
followed,  he  need  only  determine  the  correlation  between 
X*  and  r(b,  A  strong  positive  or  negative  correla¬ 

tion  will  cause  the  relative  bias  to  be  increased.  As 
we  have  discussed  earlier,  this  may  give  rise  to  lack  of 
confidence  interval  coverage. 

In  conclusion,  we  have  presented  a  two  step 
method  for  checking  variance  and  bias  increases.  Should 
the  practitioner  find  that  they  have  not  been  increased, 
he  may  feel  confident  that  the  procedure  has  worked. 

If,  however,  variance  and/or  bias  has  been  increased,  it 
is  not  clear  that  selecting  another  set  of  controls  will 
result  in  an  improvement.  In  this  case,  another  method 
should  be  considered  to  obtain  the  estimates  desired. 
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SUMMARY#  CONCLUSIONS,  AND  RECOMMENDATIONS 

This  chapter  summarizes  the  contributions  of 
this  research.  Recommendations  for  future  research  are 
also  presented. 

6.1  Research  Overview 

In  this  research  we  have  accomplished  the 
following  objectives:  (1)  to  develop  a  practical  method 
for  applying  concomitant  control  variables  in  a 
regenerative  setting  to  obtain  variance  reductions  and 
valid  confidence  intervals;  (2)  to  establish  theoretical 
properties  of  the  designated  controls;  (3)  to  validate 
and  evaluate  the  developed  method;  and  (4)  to  determine 
guidelines  for  using  the  method.  Two  types  of  standard¬ 
ized  control  variates  were  employed  in  a  two-stage 
method  to  control  both  the  numerator  and  the  denominator 
of  a  regenerative  ratio  estimator.  These  variables  have 
been  proved  to  converge  in  distribution  to  a  multi¬ 
variate  standard  normal  distribution  over  runs  (or 
cycle-batches)  of  increasing  length.  This  distribution 
served  as  a  foundation  for  performing  the  required 


regression  analyses  and  for  constructing  the  final 
confidence  interval  estimates.  In  addition,  exact 
formulas  for  the  relative  bias  and  variance  of  the  two** 
stage  estimator  were  developed. 

The  validation  and  evaluation  stage  of  this 
research  showed  that  the  developed  method  can  yield 
substantial  variance  reduction  while  maintaining  the 
confidence  interval  coverage.  Variance  reductions  of 
50%-75%  were  achieved  with  no  loss  of  coverage. 

However,  the  two-stage  method  was  observed  to  increase 
variance  under  certain  conditions. 

While  not  a  main  objective  of  this  research,  a 
new  test  for  multi-variate  normality  was  proposed.  The 
test  was  based  upon  the  univariate  Shapiro-Wilk 
statistic. 

The  results  of  this  research  confirm  the 
findings  of  others  as  to  the  practicality  of  applying 
control  variates  to  regenerative  queueing  network  simu¬ 
lations.  The  method  developed  here  would  allow  the 
practioner  to  routinely  incorporate  variance  reduction 
methods  into  his  simulations.  Following  the  guidelines 
presented  here,  some  validation  is  possible  when  the 
model  is  not  analytically  tractable. 
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The  most  important  contribution  of  this  research 
is  the  two-stage  method  and  the  expressions  for  its 
relative  bias  and  variance.  We  have  shown  these  to  be 
far  superior  to  any  other  regenerative  control  proce¬ 
dures  found  in  the  literature  today. 

6.2  Recommendations  for  Future  Research 

The  results  of  this  research  show  that  the  two- 
stage  method  is  potentially  effective  in  regenerative 
systems.  In  practice,  it  may  be  difficult  to  identify  a 
tour-defining  state  for  the  regenerative  analysis.  Even 
if  one  is  located,  the  returns  to  that  state  may  be  so 
infrequent  that  a  reasonable  number  of  tours  may  not 
occur  during  a  run  of  affordable  length.  Further 
investigation  is  needed  to  determine  how  the  two-stage 
method  performs  using  techniques  presented  in  Chapter  II 
for  approximating  a  regenerative  process. 

The  multi-variate  Shapiro-Wilk  test  presented 
needs  to  be  further  developed.  The  procedure  for 
locating  the  proper  quadratic  programming  problem  must 
be  validated.  Zn  addition,  the  null  distribution  of  the 
statistic  must  be  determined  to  obtain  tables  of  criti¬ 
cal  values.  When  fully  developed,  this  test  would  have 
widespread  applications  extending  far  beyond  simulation. 
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PROCRAM  MAIN( I NPUT, OUTPUT , TAPES* I NPUT, TAPE6-0UT PUT . TAPE?, TAPC8) 
DIMENSION  NSET ( 5000 ) 

THIS  PROGRAM  IS  DESIGNED  TO  SIMULATE  OPEN  REGENERATIVE 
QUEUEING  SYSTEMS.  TOR  EACH  CYCLE  OR  BATCH  OF  CYCLES, 

THE  TOUR  LENGTH.  TOTAL  TIME  IN  SYSTEM,  ANO  NUMBER  Or 
CUSTOMERS  AT  THE  INPUT,  BRANCH,  AND  SERVICE  POINTS  AND 
THE  ASSOCIATED  VECTORS  OF  CONTROLS  ARE  COLLECTED. 

NOTE:  TO  COLLECT  THE  BRANCH  DATA,  FOR  EACH  BRANCH  OF 
INTEREST  J.  THE  USER  MUST  COLLECT  VIA  THE  NETWORK  INPUT: 

XX(2*J-1 )»NUMBER  OF  CUSTOMERS  ARRIVING  TO  BRANCH  J 
XX(2*J) -NUMBER  OF  SUCCESSES  AT  BRANCH  J 

C0MM0N/SC0M1/ATR I B( 100), DD( 100 ) , DDL( 100) . DTNOW, I  I ,MFA,MSTOP, NCLNR 
1 , NCROR, NPRNT , NNRUN , NNSET , NTAPE . SS( 1 00 ) , SSL( 1 00 ) , TNEXT, TNOW, XX( 100) 
COMMON  QS£T(5000) 

EQUIVALENCE  ( NSET( 1 ) , QSETj 1 ) ) 

NNSET =5000 

NCR DR-5 

NPRNT-6 

NTAPE-7 

CALL  SLAM 

STOP 

END 

FUNCTION  USERF(IFN) 

C0MM0N/SC0M1/ATRIB( 100), DD( 100), DDL! 100), DTNOW, I  I ,MFA,MSTOP, NCLNR 
1 . NCROR, NPRNT, NNRUN, NNSET, NTAPE, SS( 100 ) , SSL( 100 ) , TNEXT. TNOW,XX( 100) 
C0MM0N/UC0M1 /NUMQS, CONTROL! 20 ) , NUMCUST( 20 ) , AMEAN( 20 ) , SD( 20 ) , 

1  CYSTART, SYST I ME, I  STREAM! 20 ) , NUMACT, NUMCYC, 1 8TCHSZ, I CYCCNT, NUMBR 

USERF(I)  PERFORMS  INTERARRIVAL  ON  CREATE  NODE 

USERF( 2 )  SHOULD  BE  ON  ARC  OUT  OF  CREATE  NODE  TO  CHECK  FOR  RECENER. 
USERF( 3 )  SHOULD  BE  ON  ALL  EXIT  ARCS  TO  COLLECT  TIME  IN  SYSTEM 
US£RF( J+3 ) ,  J.GT.O.  SHOULD  BE  ON  ARC  OUT  OF  QUEUE  J  TO  PERFORM 
SERVICE  TIME  ACTIVITIES 

IF  | IFN.GT.3)  CO  TO  4 
GO  TO  (1,2,3),  IFN 

INTERARRIVAL  COMPUTATIONS  (ASSUMED  EXPONENTIAL) 

ATIME-EXPON(AMEAN< 1 ), I  STREAM! 1 ) ) 

CONTROL! 1 ) -CONTROL! 1 J+ATIME 
NUMCUST ( 1 )-NUMCUST( T)+1 
USER F- AT I ME 
RETURN 

VERIFY  NEW  CYCLE  OR  NOT 

CALL  CYCLECK 

USERF-0 

RETURN 

COMPUTE  TIME  IN  SYSTEM 

SYST I ME-SYST I ME*< TNOW-ATR I B( 1 ) ) 

USERF-0 

RETURN 

PERFORM  SERVICE  TIME  COMPUTATIONS  (ASSUMEO  EXPONENTIAL) 

I- IFN- 2 

ATIME-EXPON! AMEAN( I ), I  STREAM! I ) ) 


CONTROL! I )*CONTROL( I )+ATIME 

NUMCUST ( I )*NUMCUST(  I )  +  1 

USERF*ATIME 

RETURN 

ENO 

SUBROUTINE  INTLC 

COMHON/SCOH1 /ATR I B( 1 00 ) , 00(1 00 1 , QOL< 1001. DTNOW, I  I , MFA.MSTOP, NCLNR 
1 . NCRDR, NPRNT, NNRUN, NNSET. NTAPE, SS< 100 ) . SSL( 100 > . TNEXT, TNOW.XX! 100) 
C0MM0N/UC0M1 /NUMQS, CONTROL! 20 ) f  NUMCUST( 20 ) . AMEAN(  20 ) . S0(  20 ) , 

1  CYSTART.SYSTIME, I  STREAM! 20) . NUMACT, NUHCYC, IBTCHSZ, I CYCCNT, NUMBR 
C 

C  READ  IN  NUMBERS  OF  QUEUES  AND  BRANCHES 
C 

REAO  (5.10)  NUMQS, NUMBR 
10  F0RMAT(2I2) 

NUMACT-NUMQS 

C 

C  REAO  IN  MEAN.  STANOARO  DEVIATION,  ANO  RANDOM  NUMBER  STREAM 
C  FOR  THE  INPUT  PROCESS,  SERVICE  PROCESSES,  AND  THE  BARCHINCS 
C  IN  THAT  ORDER. 

C 

00  30  1*1 . 1+NUMQS+NUMBR 

READ  (5.20)  AMEAN( I ),S0( I ), ISTREAM! I ) 

20  F0RMAT(2F10. 5,11) 

30  CONTINUE 

C  REAO  IN  BATCH  I  NO  SIZE  ANO  TOTAL  NUMBER  OF  CYCLES 
REAO  (5.35)  IBTCHSZ 
35  FORMAT  (12) 

READ  (5.37)  NUMCYC 
37  FORMAT  (16) 

1CYCCNT*0 
DO  40  1*1.20 
C0NTR0L( I )*0.0 
NUMCUST ( I )*0 
40  CONTINUE 

SYSTIME*0.0 

WR I TE( 6. 50 )  NUMQS, I BTCHSZ, NUMCYC 
50  FORMAT! IX, *THERE  ARE*, 12,*  QUEUES.  BATCH  SZ  IS  *,I2./,1X, 

1  ‘SIMULATION  IS  FOR  *, 16,*  CYCLES.*) 

WRITE<6.60) 

60  FORMAT ( IX. *FI LE*. 10X, *MEAN*. 12X. *SD*, 12X, *STREAM*) 

WRITE(6. 70)  AMEAN( 1),S0(1), I STREAM( 1) 

70  F0RMAT(1X,*INPUT*,6X, F10.5.5X, F10.5.10X, 12) 

DO  90  1*2, NUMQS+1 
J*l-1 

WRITE)  6, 80)  J.AMEAN!  I  ),S0(  I  ) ,  I  STREAM)  I  ) 

80  FORMAT! IX. ‘QUEUE  *, 12, 3X, F10. 5.5X, F10.5, 10X, 12) 

90  CONTINUE 

CYSTART-TNOW 

RETURN 

ENO 

SUBROUT  I NE  CYCLECK 

COMMON/ SC0M1 /ATR I B( 100),00( 100),0DL( 100),0TN0W, 1 1 .MFA.MSTOP, NCLNR 
1 , NCROR , NPRNT , NNRUN , NNSET , NTAPE , SS ( 1 00 ) , SSL! 1 00 ) , TNEXT, TNOW. XX( 100 ) 
COMMON/ UC0M1 /NUMQS, CONTROL! 20 ) , NUMCUST! 20 ) , AMEAN( 20 ) . S0( 20 ) . 

1  CYSTART , SYSTIME, I  STREAM! 20), NUMACT, NUMCYC, IBTCHSZ, I CYCCNT, NUMBR 
C 

C  CHECK  FOR  TOUR  DEFINING  STATE 
C 

I CHECK-0 
00  10  1*1, NUMQS 

IF  (NNQ( I )+NNACT( II  .NE.  0)  I CHECK* 1 
IF  ((CHECK  .EQ.  1)  RETURN 
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C  A  MEW  CYCLE  IS  STARTING 
C 

l*ICYCCNT 
ICYCCNT= ICYCCNT+1 
* l=MOD{ I . I BTCHS2 ) 

IF  (II  .HE.  0)  RETURN 
IF  (I  .EQ.  0)  RETURN 
C 

C  SAVE  LAST  BATCHED  DATA 
C 

ALPHA=  TNOW-CYSTART 
IF  ( IBTCHSZ  . EQ.  1 )  GO  TO  105 
DO  30  l=1,NUMQS+1 
C 

C  STANDARDIZE  SERVICE-TIME  CONTROLS 
C 

A=NUMCUST( I ) 

CONTROL( I }*(CONTROL( I )-A*AMEAN( I ) )/( S0( I )*SQRT(  A) ) 

30  CONT I NUE 

DO  40  1=1 , NUMBR 
J=  I +NUMQS-M 
NUMCUST ( J ) *XX(  I *2-1 ) 

C 

C  STANDARDIZE  FLOW  CONTROLS 
C 

A=XX( 2*1-1 ) 

CONTROL! J )=( XX( 2* I )-A*AMEAN( J ) )/( SD( J )*SQRT( A) ) 

40  CONT I NUE 

105  WRITE! 8,107)  ALPHA, SYST IME, ( NUMCUST! I ) .CONTROL! I ), 1=1 ,4) 
107  FORMAT! 2F1 0.4, 4 ( 1 6, F10. 4 ) ) 

C 

C  CHECK  FOR  END  OF  SIMULATION 
C 

I F( ICYCCNT  .G£.  NUMCYC)  MST0P=-1 
C 

C  REINITIALIZE  FOR  NEXT  CYCLE 
C 

SYST I ME=0 . 0 

DO  120  1=1, NUMQS* 1 +NUMBR 
CONTROL! I )=0.0 
NUMCUST! I )=0 
120  CONTINUE 

CYSTART»TNOW 

RETURN 

END 


CEN, DENNY, MM1  QUEUE, 6/15/81 ; 
UM.1,2.500; 

NETWORK; 

CREATE, USERF(  1 ) , ,  1 ; 
ACT,  USERF(  2 ) ; 

QUEUE! 1 ) ; 

ACT/ 1 , USERF  ( 4 ) ; 

COON; 

ACT, USERF(  3); 

TERM; 

ENONETWORK; 

INIT.O; 

FIN; 

1  0 

1,0  1.0  9 

.50  .50  9 

15 

375001 
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PROGRAM  MAIN< I NPUT, OUTPUT, TAPE5= INPUT, TAPE6*0UTPUT, TAPE7, TAPES ) 

THIS  PROGRAM  IS  DESIGNED  TO  SIMULATE  CLOSED  RECRENERAT I VE 
QUEUEING  SYSTEMS.  FOR  EACH  CYCLE  AR  BATCH  OF  CYCLES,  IT 
COLLECTS: 

1 )  TOUR  LENGTH  ALPHA 

2)  TOTAL  CUSTOMER  TIME  IN  SYSTEM 

3)  NUMBER  OF  COMPLETED  CUSTOMERS 

4)  BUSY  TIME  FOR  EACH  SERVICE  ACTIVITY 

5 )  NUMBER  OF  CUSTOMES  ARRIVING  TO  EACH  BRANCH  OR  ACT. 

6)  CONTOL  VARIABLE  FOR  EACH  BRANCH  AND  ACTIVITY 

NOTE:  IN  THE  NEWORK,  THE  USER  MUST  COLLECT 

XX( 2*J- 1 )=THE  NUMBER  OF  CUSTOMERS  ARRIVING  TO  BRANCH 
POINT  J 

XX( 2*J )=THE  NUMBER  OF  SUCCESSES  AT  BRANCH  J 
DIMENSION  NS£T( 5000 ) 

COMMON/SCOM1/ATR I B< 100 ) , DD( 100 ) , 0DL( 100 ) , DTNOW, I  I , MFA, MSTOP, NCLNR 
1 , NCRDR, NPRNT . NNRUN, NNSET, NTAPE, SS( 100 ) , SSL( 100 ) , TNEXT, TNOW,XX( 100 ) 
COMMON  QSET ( 5000 ) 

EQUIVALENCE  (NSET( 1 ),QSET( 1 ) ) 

NNSET*5000 

NCR0R=5 

NPRNT*6 

NTAPE=7 

CALL  SLAM 

STOP 

END 

FUNCTION  USERF(IFN) 

COMMON/SCOM 1 /ATR I 8( 1 00 ) . 00( 1 00 ) , DOL( 1 00 ) , DTNOW, I  I , MFA, MSTOP, NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE . SS ( 1 00 ) , SSL( 1 00 ) , TNEXT , TNOW, XX( 1 00 ) 
COMMON/UCOM 1 / NUMQS , CONTROL ( 20 ) , NUMCUST( 20 ) , AMEAN( 20 ) , SD( 20 ) , 
1CYSTART, SYST I  ME, I STREAM( 20 ), NUMBR, NUMCYC, IBTCHSZ, ICYCCNT, NUMCSTS, 

1 NUMCOMP 

IN  NETWORK,  USERF( I  )  IS  ON  ARC  LEAD  I NC  INTO  QUEUE  I. 

USERF( l+NUMQS)  IS  ON  ARC  OUT  OF  QUEUE  I. 

IF  (IFN  .GT.  NUMQS)  CO  TO  4 

VERIFY  NEW  CYCLE  OR  NOT  ANO  COMPUTE  TRAVERSE  TIME 

IF  (IFN  .NE.  1)  RETURN 

SYST  I  ME=*SYST  I  ME>(  TNOW-ATR I  B(  1  )  ) 

NUMC0MP*NUMC0MP+1 
CALL  CYCLECK( IFN) 

USERF-0 

RETURN 

PERFORM  SERVICE  TIME  COMPUTATIONS 


l»  I FN-NUMQS 

AT IME*£XPON( AMEAN( I ), ISTREAM( I ) ) 

CONTROU I )*CONTROL( I J+ATIME 

NUMCUST( I )*NUMCUST(  I )  +  1 

USER F* AT  I ME 

IF  ( I  .NE.  1 )  RETURN 

ATRIB(1)*TN0W+ATIME 

RETURN 

ENO 

SUBROUTINE  INTLC 

C0MM0N/SC0M1/ATRI8< 100),00( 100),00L( 100),0TN0W, » I , MFA, MSTOP, NCLNR 
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1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE , SS( 1 00 ) . SSL ( 1 00 ) . TNEXT . TNOW, XX( 100 ) 
COMMON/UCOM1 /NUMQS , CONTROL( 20 ) , NUMCUST( 20 ) . AMEAN( 20 ) .  SD<  20 ) . 

1 CYST ART, SYSTIME, I STREAM! 20 ) , NUMBR, NUMCYC, IBTCHSZ, ICYCCNT, NUMCSTS, 

1 NUMCOHP 
C 

C  READ  NUMBERS  OF  QUEUES.  BRANCHES.  ANO  CUSTOMERS 

A  C 

READ  (5,10)  NUMQS.NUMBR, NUMCSTS 
10  FORMAT (314) 

C 

C  READ  MEANS,  STANOARO  DEVIATIONS,  ANO  RANDOM  NUMBER  STREAMS 
C  FOR  QUEUES  THEN  BRANCHES 

,  C 

*  DO  30  1=1, NUMQS+NUM8R 

READ  (5,20)  AMEAN( I ) , SD( I ) , I  STREAM! I ) 

20  F0RMAT(2F10.5, 11 ) 

30  CONTINUE 

C 

C  READ  BATCH  SIZE  AND  TOTAL  NUMBER  OF  CYCLES 
C 

READ  (5,35)  IBTCHSZ 
35  FORMAT  (12) 

READ  (5,37)  NUMCYC 
37  FORMAT  ( 16) 

ICYCCNT=0 
DO  40  1=1.20 
CONTROL! I }=0.0 
NUMCUST( I )=0 
40  CONTINUE 

SYSTIME=0.0 
NUMCOMP-O 
DO  45  1=1,2 
45  XX(l)=0.0 

WR I TE( 6, 50 )  NUMQS, I BTCHSZ, NUMCYC, NUMCSTS 
50  FORMAT! IX, *THERE  ARE*. 12,*  QUEUES.  BATCH  SZ  IS  *, 12,/, IX, 

1  *SIMULATION  IS  FOR  *, 16,*  CYCLES. *,/, IX, 

1  *THERE  ARE*, 13,*  CUSTOMERS.*) 

WRITE(6,60) 

60  FORMAT! 1 X. *F I LE*. 10X, *MEAN*, 1 2X, *S0*, 1 2X, “STREAM* ) 

00  90  1=1, NUMQS 

WR  I  TE(6, 60)  I.AMEA N(  I  ),S0(  I  ),  I  STREAM(  I  ) 

80  F0RMAT(1X,*QUEUE  *, 12, 3X, F10. 5, 5X, F10. 5, 10X, 12) 

90  CONTINUE  * 

CYSTART=TNOW 

C 

C  LOAO  INITIAL  CUSTOMERS 
C 

ATRIB{ 1 )*0.0 
00  100  1=1, NUMCSTS 
100  CALL  FILEM(1,ATRIB) 

RETURN 

ENO 

SUBROUTINE  CYCLECK(IFN) 

COMMON/SCOM1 /ATR I B( 1 00 ) , D0( 1 00 ) , DOL( 100 ) , OTNOW, I  I , MFA, MSTOP, NCLNR 
«  1 , NCROR .NPRNT, NNRUN , NNSET , NTAPE , SS( 1 00 ) , SSL ( 1 00 ) , TNEXT, TNOW. XX( 1 00 ) 

COMMON/UCOMI /NUMQS, CONTROL! 20 ) , NUMCUST( 20 ) , AMEAN( 20 ) , S0( 20 ) , 
1CYSTART, SYSTIME, I  STREAM! 20), NUHBR, NUMCYC, IBTCHSZ, ICYCCNT, NUMCSTS, 

1 NUMCOHP 

DIMENSION  BUSY( 20) 

C 

„  C  CHECK  FOR  TOUR  DEFINING  STATE 

•  C 


l*NNQ( 1 )*NNACT( 1 )♦! 
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IF  (I  .NE.  NUMCSTS)  RETURN 
A  NEW  CYCLE  IS  STARTING 


C 

C 

C 


I = I CYCCNT 
I CYCCNT* ICYCCNT+1 
I  I *MOD( I, IBTCHSZ) 

IF  (II  .NE.  0)  RETURN 

IF  (I.EQ.O  .AND.  IBTCHSZ.NE. 1 )  RETURN 


SAVE  LAST  BATCHEO  DATA 


20 


C 

C 

C 


ALPHA=TNOW-CYSTART 
00  20  1*1, NUMQS 
CALL  TTUTL( I , UAVG, U I  NT, OT ) 
BUSY( I )=UINT 
CONTINUE 

IF  (IBTCHSZ  .EQ.  1)  GO  TO  105 
DO  30  1*1, NUMQS 


STANDARDIZE  SERVICE-TIME  CONTROLS 


30 


CONTROL ( I )=( CONTROL! I )-A*AMEAN( I ) )/( SD( I )*SQRT( A) ) 
CONTINUE 
DO  40  1*1, NUMBR 
J*1+NUMQS+I 
NUMCUST ( J )=XX( 1*2-1) 


C  STANDARDIZE  FLOW  CONTROLS 
C* 

CONTROL! J)=(XX( I *2 ) -A*AM£AN( J ) )/(SD(J )*SQRT(A) ) 
BUSY(J)=0.0 

WR"TE?»fl07)  ALPHA, SYST IME, NUMCOMP, 

1  ( BUSY( I ) , NUMCUST( I ) , CONTROL! I ), 1*1,6) 
F0RMAT(2F10.4, I4,6(F10.4, |4,F10.4)) 


40 

105 


107 
C 
C 
C 


CHECK  FOR  END  OF  SIMULATION 

I F ( I CYCCNT  .GE.  NUMCYC)  MSTOP-- 1 


C  REINITIALIZE  FOR  NEXT  CYCLE 
C 

SYST I ME*0 . 0 
NUMC0MP*0 

DO  120  1*1, NUMQS+NUMBR 
CONTROL! I )*0.0 
NUMCUST ( I )»0 
CONTINUE 
00  125  1*1.2 
XX( I )*0.0 
CYSTART*TN0W 
CALL  CLEAR 
RETURN 
ENO 

SUBROUT  I NE  TTUTL( I ACT , UAVG, U I  NT, DT ) 

DIMENSION  NSET ( 1 ) 

DIMENSION  NN0(2,2) 

COMMON / SC 0M1  / * ATR I B(  100) , D0(  100 ) , DDL(  100 ), OTNOW,  '  ' 

1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE, SS( 1 00 ) , SSL( 100), TNEXT , TNOW, XX(  1 00 ) 


120 

125 


H  I'm  'iV  TMiinfr&ii 
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IF  (NSET(K)-II)  90,80,180 
80  NFSN*NSET( K+6 ) 

I F( NFSN.GT.O.OR. NSET( K+7 ) . CT.O )  CO  TO  180 
NTYPal 
KK*K+7 
GO  TO  100 
90  NSSR»NSET ( K+4  J 

IF  (NSSR. EQ.O)  GO  TO  180 
NTYP*2 

KK»K+-6+NSET ( K+5 ) 

C 

C***  COMPUTE  ANO  WRITE  SERVER  STATISTICS 
C 

100  IF  (IK.EQ. 1)  GO  TO  110 
I  K*1 

110  NCAP*NSET ( KK+1 3 ) 

I NOEX=NSET ( KK+3 ) 

IF  ( INDEX. NE. I  ACT)  GO  TO  170 
IF  (XT.LE.O. )  GO  TO  200 
XBUSY=NSET( KK+12 ) 

I BUSY=XBUSY 
BLCK3NSET< KK+1 1 ) 

TDEL3TNOW-QSET(  KKffi ) 

XBT3XBUSY*TDEL 
U I NT*(QSET( KK+4 ) +XBT ) 

UTIL3  U I  NT/XT 
UAVG  3  UTIL 
OT  «  XT 

STD3 ( QSET( KK+14) +XBUSY*XBT )/XT-UT I L*UT I L 
STD-S I GN( SQRT( ABS( STD )) , STO ) 

BLCK3(QSET( KK+5)+BLCK*T0EL)/XT 
XBMAX=OSET( KK+7 ) 

X I MAX=QSET ( KK+8 ) 

IF  ( NCAP.GT. 1 )  GO  TO  130 
IF  ( I  BUSY , EQ . 1 )  GO  TO  120 
IF  (TDEL.GT.XIMAX)  XIMAX=TDEL 
GO  TO  140 

120  IF  ( TDEL.GT.XBMAX)  XBMAX=*TDEL 
GO  TO  140 

130  IF  (XBUSY.CT.XBMAX)  XBMAX3XBUSY 
XI0L*-NSET(KK) 

IF  (XIDL.GT.XIMAX)  XIMAX-XIOL 
140  CONTINUE 
GO  TO  200 
170  KK3NSET(KM-1 ) 

IF  (KK.GT.O)  GO  TO  100 
180  CONTINUE 
GO  TO  200 

190  CALL  ERROR  (568) 

200  RETURN 


SUBROUTINE  UMONT( ITRACE) 

COMMON/SC 0M1/ATR I B( 1 00 ) , 00( 1 00 ) , 00L( 1 00 ) , OTNOW, I  I , MFA, MSTOP, NCLNR 
1 , NCROR, NPRNT, NNRUN, NNSET, NTAPE, S$( 100 ) , SSL( 100 ) , TNEXT , TNOW, XX( 100 ) 
XX< 50)*NNACT( 1 )+NNQ(  1 ) 

XX(51 )3NNACT(2)+NNQ(2) 

XXI 52 )»NNACT( 3 )*NNQ( 3 ) 

RETURN 

ENO 
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ft 


* 


CEN.OENNY. SYSTEM  4,6/15/81; 
LIM. 4,2, 10; 

NETWORK; 

ONE  COON; 

ACT,USERF(1 ); 

QUEUE( i ); 

ACT(5)/1,USERF(5); 
ASSIGN,XX( 1 )=XX( 1 )+1 . ; 
ACT, ,.25,02; 

ACT, , . 75.Q3A; 

Q2  QUEUE(2); 

ACT/2,  USERF(  6  )..Q4; 
Q3A  ASSIGN, XX(2)=*XX(2)-M.; 

ACT,,,Q3; 

Q3  QUEUE!  3); 

ACT/3.USERF! 7) , ,Q4; 

Q4  QUEUE! 4 ) ; 

ACT/4, USERF( 8); 
ASSIGN, XX! 3 )*XX( 3 )+1 . 0* 
ACT,,, 9, ONE; 

ACT, , . 1 ,Q3B; 

Q3B  ASSIGN, XX{4)-XX!4)+1; 

ACT,,,Q3; 

ENDNETWORK; 
lit  IT,  0; 


FIN; 

4 

4  7 

10. 

10. 

9 

1.5 

1.5 

9 

1.0 

1.0 

9 

.5 

.5 

9 

.75 

.433 

9 

.1 

.3 

9 

24 

60000 


E 


■* 
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» 


5 

10 

20 

C 

C 

C 


C 

c 

c 


25 


30 

C 

C 

C 


40 


50 


60 

70 


C 

C 


PROGRAM  REGRE( INPUT, OUTPUT. TAPE5* INPUT. TAPE6*0UTPUT, TAPE?) 
DIMENSION  X( 50 ) ,CB( 50, 1 ) . Y(50  ),CT(50,4 ) ,Z( 50) ,CB8AR( 1 ), 

1  CTBAR(4),SICMAB(1,1 ) ,  SIGMAT(4, 4 )  ,RR(  50,2 ), 

1  DELTAS! 1 ) , DELTAT! 4 ) . CAMMA( 1 ) , BETA! 4 ) . 

1  A(50 ) ,8( 50 ) , S IGMABI (1,1) , SIGMAT 1(4,4) 

NTC0NT*4 

NBCONT* 1 

I SAMPSZ*50 

ZVAL* 1.645 

NEXPCT-0 

NUMEXP*1 

RTRUE*1.0 

B I ASR*0 . 0 

BIASRK-0.0 

BIASRIIT-0.0 

8 I ASRHB*0 . 0 

VARR=0. 0 

VARRHAT*0.0 

DO  20  1*1, ISAMPSZ 

REA0(5, 10)  X( l),Y( l),CB{ l,1).CT( 1,1) 

F0RMAT(2F10.4,2( 10X, F10.4) ) 

CONTINUE 

PERFORM  BOTTOM  REGRESSION 

CALL  MEAN(XBAR,CBBAR, I SAMPSZ , NBCONT, X, CB) 

CALL  SIGMA(C8,C8BAR, ISAMPSZ, NBCONT. SIGMAS) 

CALL  OELTA( X, CB , XBAR , CBBAR, I SAMPSZ. NBCONT, DELTAS) 

CALL  RECRESS( GAMMA, S I GHA8, S I GMAB I , DELTAB, NBCONT ) 

DETERMINE  DENOMINATOR 

BBAR*0.0 

DO  30  1*1. ISAMPSZ 
SUM*0 . 0 

DO  25  J*1, NBCONT 

SUM*SUM+GAMMA( J )*C8( I , J ) 

B(  I  )*X(  I  J-SUM 
BBAR=BBAR+B( I ) 

CONTINUE 

8BAR-B8AR/ I SAMPSZ 
PERFORM  TOP  REGRESSION 


YBAR*0 . 0 

DO  40  1-1, ISAMPSZ 
YBAR*YBAR+Y( I ) 

Y8AR*YBAR/ I SAMPSZ 
R-Y8AR/XBAR 
ZBAR*Y8AR~R*B8AR 
00  50  1*1, ISAMPSZ 
Z<  I  J«Y<  I  )-R*B(  I ) 

DO  70  1*1, NTCONT 
CTBAR( I )*0.0 
DO  60  J*1. ISAMPSZ 

CTBAR( I )*CTBAR( I )*CT( J, I ) 

CTBAR( I )*CT8AR( I )/ I SAMPSZ 
CONTINUE 

CALL  S I OMA(CT,CTBAR, I SAMPSZ, NTCONT, SIGMAT) 

CALL  DELTA(Z,CT,ZBAR,CTBAR, I SAMPSZ. NTCONT, DELTAT) 
CALL  RCCRESS( BETA, S I GMAT , S I  GHAT  I , DELTAT , NTCONT ) 


DETERMINE  NUMERATOR 


ono  ooo 
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c 


75 

80 


ABAR-0.0 

DO  80  l“1 » ISAMPSZ 
SUM-0.0 

00  75  J-1.NTC0NT 

SUM*SUM+BETA< J )*CT( I , 
A( I )»Y< I )-SUM 
A8AR—ABAR+ A ( I ) 

CONTINUE 

ABAR-ABAR/ ISAMPSZ 


J) 


GET  POINT  ESTIMATE 


RHAT-ABAR/BBAR 
RHATB* YBAR/ 8BAR 
RHATT-ABAR/XBAR 

GET  CONF I DENCE  INTERVAL 


90 

100 


110 


120 


SUM1-0.0 

SUM2-0.0 

DO  90  1-1 , NTCONT 

SUM1*SUN1+BETA( I )*CTBAR(I ) 

Dy  100  I  —  1 , NBCONT 

SUM2*SUM2+CAMHA( I )*C88AR( I ) 
SUMTRM-SUM1 -RHAT*SUM2 
VAR* 0.0 
VAR 1-0.0 

°°VIr2vArJ{ A?n^RMAT*B{  I  )*SUMTRM )**2 
VAR1-VARJ*!  Y(  I  )-R*X|  » ) )**2 
CONTINUE 
AMf»  I  SAMPSZ 

VAR-VAR/<AN*BBAR**2*( AN- 1.0)) 
VAR1-VAR1 / ( AN*X8AR**2*( AN-1  - 0 ) ) 

2VAL 1 -ZVAL*SQRT ( VAR) 

2'VAL2»ZVAL*SQRTi  VAR1 ) 

CIL0W-RHAT-2VAL1 

ClMIGH-RHAT+ZVALI 

CIL0W1-R-ZVAL2 

CIHI0H1-R+ZVAL2 

br-r-rtrue 

BRHAT-RHAT-RTRUE 

BRHATB-RHATB-RTRUE 

SWIDTrn'Kctow.oi.^.v-,... 

FORMAT! 1X.8F8. 5) 

Bl ASR-BI ASR+BR _ 

8 I ASRH-B I ASRH+BRHAT _ 

B I ASRHB-B I ASRHB+BRHATB 
8 I ASRNT-B I A5RHT+BRHATT 
VARR- VARR+VAR 1 
VARRNAT-VARRHAT+VAR 

S»S5S!^. 

TCffiSIMfcUw «.  ™ » 

an-numcxp 

81 ASR-BI A8R/AN 
BIA8RN»8IASRH/AN 
biasrmt-biasrht/an 
8 1  A8RH8-B 1 A8AKB/AN 
RH1-0.0 


CIL0W1.CIHICH1.VAR1 
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140 


10 

20 


20 

10 


10 

20 

30 


10 


20 

30 

40 


RH2-0.0 
VAR1-VARR/AN 
VAR2-VARRHAT/AN 
VARRED- ( VAR 1 -VAR2 ) /VAR 1 

WRITE(7, 140)  B I ASR , B I ASRH , B I ASRHB , B I ASRHT , VARRED 
FORMAT ( IX, ‘BIAS  IN  CLASSICAL-*, F10. 4./, 

1  1X,*BIAS  IN  CONTROLLED  TOP  AND  BOTTOM-*, FI 0.4, /, 

1  IX, *BI AS  IN  CONTROLLED  BOTTOM-*, F10. 4,/, 

1  1X,*BIAS  IN  CONTROLLED  TOP-*, F10.4,/, 

1  IX, *VAR  REDUCTION  IN  CONTROLLED  TOP  ANO  BOTTOM*. FI 0.4) 


STOP 

ENO 

SUBROUT  I NE  DELTA( X, C, XBAR . CBAR , N, NP, DELTAMT ) 

DIMENSION  X(N),C(N,NP).CBAR(NP),DELTAMT(NP) 

AN-N-1 

DO  20  1-1, NP 
SUM-0.0 
00  10  J*1 ,N 

SUH-SUM+(C(J, I )-CBAR( I ) )*(X( J  J-XBAR) 

DELTAMT ( I ) -SUM/ AN 
CONTINUE 
RETURN 
ENO 

SUBROUT  I NE  RECRESS( BETA, S I CMAMT, S IGMAI , DELTAMT, NP) 

DIMENSION  BETA(NP),SIGMAMT(NP,NP),D£LTAMT(NP),SIGMAI(NP,NP) 
CALL  MAT INV( SI CMAMT, NP, SIGMA I ) 

DO  10  1-1, NP 
BETA( I )-0.0 
DO  20  J-1.NP 

8ETA( I )-BETA< I J+SICMAI { I , J )*DELTAMT( J ) 

CONTINUE 

RETURN 

END 

SUBROUTINE  SIGMA(C,C8AR,N,NP,SICMAMT) 

DIMENSION  C( N, NP ) , CBAR( NP) , SIGMAMTJ  NP.NP) 

AN-N-1 

DO  30  1-1, NP 
DO  20  J-1.NP 
SUN-0.0 
DO  10  K-1.N 

SUM-SUM+(C(K,  I  )-CBAR(  I  )  )*(  C(  K,  J  )-CBAR(  J ) ) 

SI CMAMT ( I , J )-SUM/AN 
CONTINUE 
CONTINUE 
RETURN 
ENO 

SUBROUTINE  MEAN( XBAR, CBAR, N, NP,X,C) 

DIMENSION  CBAR( NP),X(N),C(N,NP) 

AN-N 


XBAR-0.0 
DO  10  1-1, NP 
C8AR( l)-0.0 
00  30  J-1.N 
XBAR-X8AR+X( J ) 

DO  20  1-1, NP 

C8AR( I )>CBAR( I )*C( J, I ) 
CONTINUE 
X8AR-XBAR/AN 
00  40  1-1, NP 

C8AR( I )-CBAR( I )/AN 
RETURN 
ENO 

SUBROUTINE  MATINV(AMT,N,AINV) 


0 
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